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ABSTRACT 


kl 

A  computer  program  is  described  for  use  on  the  IBM-704/7090 
electronic  data  processing  machine  or  any  large  computer  accepting 
FORTRAN.  The  necessary  modifications  for  use  on  these  two  com¬ 
puters  are  simple  and  self-evident.  The  computer  program  permits 
the  computation  of  detailed  ray  patterns  portraying  the  spatial  dis¬ 
tribution  of  rays  emitted  from  a  transmitter  whose  geographic  co¬ 
ordinates  with  respect  to  the  center  of  the  earth  are  known.  This 
program  deals  with  the  solution  of  the  differential  equations,  given 
by  Hamiltonian  optics,  for  ray  paths  in  non-isotropic,  three- 
dimenaionally  nonhomogeneous  media  whose  complex  phase  refractive 
index  is  given  by  the  Appleton-Hartree  formula. 


in  presenting  ,$n 

account  of  the  current  status  of  the  development  of  this  program, 
which  has  yielded  many  useful  results,*  Pre  sente  del  so  are^gample 
calculations  as  well  as  some  results  that  hftVe  isisn  obtained  by  using 
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SECTION  I 

INTRODUCTION 


What  is  customarily  referred  to  as  reflection  in  electromagnetic  wave 
propagation,  is  actually  the  result  of  an  integrated  effect  of  a  pheno¬ 
menon  of  refraction  (i.  e. ,  the  gradual  changing  of  the  direction  of 
the  electromagnetic  energy  transport  vector).  Insofar  as  the  pheno¬ 
menon  of  refraction  is  concerned  it  is  well  established  that  the 
spatial  gradient  of  electron  density  plays  a  crucial  role  in  control¬ 
ling  the  propagation  behavior  of  an  electromagnetic  wave.  Although 
this  is  well  known,  propagation  studies  so  frequently  incorporate  a 
curious  mixture  of  simple  refraction  phenomena  (in  the  use  of  Snell's 
law)  and  of  mirror-like  reflection.  Hence,  in  these  studies  they 
omit  an  accounting  of  these  spatial  gradients  of  electron  density  as 
well  as  their  variation.  As  generally  known,  these  oversimplified 
studies  derive  from  Snell's  law  the  condition  for  reflection,  that 
being,  that  at  the  spatial  point  of  reflection,  the  electron  density  must 
attain  the  value  of  1.  24  x  10*  f 2  cos^  I  electron/ cc,  where  f  is  the 
propagation  frequency  in  Me/ sec  and  I  is  the  angle  of  incidence  of 
the  ray  upon  the  first  layer  of  the  assumed  stratified  ionosphere. 

Such  a  simplified  approach  is  necessary  because  of  the  difficult  task 
that  leads  to  a  numerical  solution  of  a  propagation  problem  which 
incorporates  these  electron  density  gradients  in  its  solution. 

With  increasing  use  of  transmitters  in  satellites,  as  well  as,  for 
the  understanding  of  the  behavior  of  wave  propagation  under  severly 
abnormal  atmospheric  conditions,  it  becomes  important  to  take  a 
realistic  account  of  the  spatial  gradient  of  the  electron  density  in 
electromagnetic  wave  propagation  studies.  Only  in  this  manner  will 
it  become  possible  to  usefully  utilise  the  new  satellite  propagation 
techniques  in  studies  designed  toward  the  understanding  of  the  atmos¬ 
pheric  ionisation-deionisation  phenomena  and  through  this  the  detailed 
structure  of  the  ionosphere. 

Within  the  last  ten  years  some  effort  has  been  made  in  constructing 
analogue  computers  for  the  study  of  ray  propagation  which  accounted 
for  spatial  electron  density  gradients  (as  for  example  D.  F.  Hartree, 
et  al. ,  Manchester,  England;  M.  S.  Wong,  AFCRL,  Bedford,  Mass. ). 
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Some  of  these  analogue  computers  were  and  still  are  limited  to  spatial 
electron  density  gradients  in  a  particular  direction  thus  forcing  the 
propagation  problem  into  a  two-dimensional  plane,  or  to  the  study  of 
the  behavior  of  refraction  on  wave  propagation.  These  constraints 
are  built  into  the  analogue  computer  and  are  not  easily  changed. 

One  approach  that  avoids  these  constraints  is  the  use  of  a  large 
electronic  data  processing  system  where  the  ordered  logical  flow 
controlling  any  calculation,  is  easily  varied.  Combining  such  a 
programmed  computer  with  Hamiltonian  optics,  which  give  the  desired 
ray  tracing  equations  for  a  nonhomogeneous ,  non- isotropic  mediu, 
and  the  Appleton- Hartree  formula  for  the  complex  refractive  indes, 
permits  in  addition  to  three  dimensional  ray  tracing,  the  simultaneous 
study  of  numerous  other  variables  of  the  propagation  problem. 

Such  an  approach  to  the  ray  trace  propagation  problem  is  presented 
in  its  present  state  of  development.  A  great  deal  of  improvement  in 
some  of  the  routines  is  possible.  As  a  result,  the  writer  would  like 
to  encourage  correspondence  concerning  these  matters.  In  addition, 
it  should  be  stated,  that  a  computer  program  for  solving  the  three- 
dimensional  ray  trace  problem  has  also  been  written  for  the  Ferranti 
Mercury  Computer  at  Manchester  University  by  C.  B.  Haselgrove  and 
J.  Haselgrove. ^ 
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SECTION  II 

COMPUTATIONAL  PROCEDURE 


A,  RAY  TRACE  EQUATIONS 


When  electromagnetic  waves  are  propagated  through  a  medium 
whose  dielectric  constant  or  index  of  refraction  is  a  varying  function 
of  the  path,  the  waves  undergo  a  change  in  direction,  or  refractive 
bending,  and  a  retardation  in  the  velocity  of  propagation.  The  spatial 
relationship  expressing  this  angular  bending  of  a  ray  of  an  electro¬ 
magnetic  wave  can  be  determined  by  basing  the  theory  of  rays  and 
waves  on  a  variational  principle  (Fermat's)  in  space.  By  a  ray  is 
meant  the  path  travelled  by  the  transport  vector  of  electromagnetic 
wave  energy  between  the  transmitter  and  an  associated  electromagnetic 
field- intensity  point  in  space.  This  Hamilton  Theory  starts  from 
the  variational  principle  b  Jpds  =  0,  where  P  is  a  medium  function 
or  index  of  refraction,  depending  on  position  and  direction.  From 
this  principle  the  theory  constructs  the  properties  of  systems  of 
rays  and  the  waves  associated  with  them  (i.  e. ,  extremals  and  trans¬ 
versals,  in  the  language  of  the  calculus  of  variations).  Because  of 
the  stationarity  in  time,  the  theory  may  be  regarded  as  a  statical 
one,  the  rays  being  fixed  curves  in  space  and  the  waves  fixed  surfaces. 
Neither  wave-length  nor  frequency  is  involved.  Likewise  the  waves 
form  a  continuous  set  of  surfaces,  not  distinguished  as  crests  and 
troughs.  This  theory,  whether  in  the  form  preferred  by  Hamilton 
or  otherwise,  has  been  the  subject  of  many  books  under  the  general 
title  "Geometric  Optics". 


Thus,  applying  Hamiltonian  optics  leads  to  the  general  Hamilton's 
Equations^  for  ray  paths  of  electromagnetic  waves  in  a  three- 
dimensional  non-isotropic  nonhomogeneous  medium.  From  them 
Haselgrove^  has  derived  the  following  set  of  equations  for  ray  paths 
in  a  spherical  coordinate  system  in  a  format  suitable  for  numerical 
integration  on  high  speed  computers: 


dr 

dt 


(1) 
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(6) 


In  these  equations  r,  0,  and  9  are  the  spatial  coordinates  of  a 
spherical  system;  P  is  the  arbitrary  index  of  refraction;  3  is  a  vector 
directed  normal  to  the  phase  fronts  of  the  ray  of  magnitude  P  with  dr 
9,  and  0^  its  respective  components  in  the  r,  6,  and  9  directions; 

X  is  the  time  of  phase  travel  along  the  ray  path  (i.  e.  ,  f  Ax/c  =  the 
number  of  wavelengths  in  the  medium  along  the  ray  path,  where  f 
is  the  electromagnetic  wave  frequency  and  c  the  velocity  of  light 
in  vacuum). 

It  is  noteworthy  that  the  partial  derivatives  of  P  appear  explicitly  in 
Equations  1  to  6,  in  accordance  with  the  fact  that  the  gradients  of  p 
play  crucial  roles  in  determining  the  spatial  ray  paths. 

This  closed  set  of  first-order  partial  differential  equations  which 
will  describe  the  propagation  behavior  of  an  electromagnetic  wave 
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under  geometric  optics  conditions,  can  be  integrated  simultaneously 
by  a  point-by-point  numerical  process  if  expressions  can  be  develop¬ 
ed  for  the  necessary  derivatives  of  the  phase  refractive  index  P. 

The  quantity  P  and  its  derivatives  are  obtained  by  using  the  Appleton- 
Hartree  formula 3  as  the  definition  of  the  complex  phase  refractive 
index  M. 

The  derivatives  are  derived  under  the  conditions  of  ray  optics,  that 
is,  that  the  imaginary  part  of  M2  is  very  much  smaller  than  the  real 
part.  As  an  aid  for  computer  use  and  comparison  with  published 
works  of  others^.^  the  Appleton-Hartree  formula  is  written  as: 

m2  =  i,  -  j«>2  =  i  -  (7| 


D  =  2(1  -  X  -  jZ)  <  1  -  jZ)  -  Y2  sin2  *  +  S  (8) 

S  =  ±[  (Y  sin*)4  +  4Y2(1  -  X  -  jZ)2  cos2*  }  (9) 

where 

M  =  the  complex  phase  refractive  index 

X  =  a  scalar  quantity  representing  the  normalized  electron 
density 

2  oo  ^ 

4irNe  p 

2  =  2 
moo  oo 

oo^  =  plasma  frequency  at  the  spatial  point 
oo  =  angular  wave  frequency  =  2irf 
m,  e  =  mass  and  charge  of  an  electron 
N  =  electron  density  at  a  spatial  point 
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Y  =  normalized  magnitude  of  the  earth's  magnetic  field 
vector  ?  = 


mcCO  to 


05  =  magnitude  of  the  gyromagnetic  frequency  of  the  electron 

about  the  earth's  magnetic  field 

Z  =  normalized  collision  frequency  =  (V/05) 

V  =  collision  frequency  at  a  spatial  point  in  collisions  per 
second 

=  angle  defined  by  the  inner  product  of  the  vectors  and 
?  = 


cos 


•*[■ 


Vr  ♦  %Ye  +_W_  j 


(HY) 


ck 


K  =  —  =  imaginary  part  of  the  complex  phase  refractive 
index 


c  =  velocity  of  light  in  vacuum 

k  =  absorption  coefficient  of  the  wave  per  unit  length  of 

path  (it  is  proportional  to  the  conductivity  of  the  medium) 


It  is  noted  that  there  are  two  possible  values  for  the  complex  index 
of  refraction  M  corresponding  to  the  plus  and  minus  sign  on  S 
which  represent  two  different  modes  of  ionospheric  propagation. 

These  are  commonly  called  "ordinary"  and 'fextraordinary"modes 
for  the  plus  and  minus  sign  respectively.  Also,  the  Appleton-Hartree 
formula  (Equations  7  to  9)  is  notable  in  that  it  presents  p,  which 
actually  is  a  spatial  function  of  six  variables,  in  the  form  containing 
purely  algebraic  operations  on  factors  or  terms  each  of  which  is  a 
function  of  at  most  three  variables,  that  is,  either  of  r,  9,  9  or  0  , 
Og,  0^  .  This  reduces  the  representation  of  P  to  a  numerical 
problem,  easily  solvable  to  current  computer  techniques. 
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If  i  represents  any  one  of  the  spatial  spherical  coordinates  r,  0, 
and  t,  then  the  partial  derivatives  of  n  with  respect  to  the  com¬ 
ponents  of  the  wave  normal,  $  ,  can  be  easily  shown  to  be 


a  a  a j,  -  /  0.  Y  COS  t  -  0  Y. 

5|i  _  3n  8t  _  an  j  1 _ _ _ i_ 

do.  =  a^  ao.  =  a^M  .  * 

i  i  \  0  Y  sint 


(10) 


This  useful  transformation  also  enjoys  the  following  property: 


w  wr- 

When  t  -*  o,  -»  o,  ^ — >  oo  but  o 


at 


an 


80. 


30. 


(ID 


To  evaluate  8n/  8o^;  the  necessary  partial  derivative  is: 


a?  _  R  8M  _ 

0^  “  Re  Tf  "  R 


(M2-!)  (Y2sintcost) 
MD 


.  an 
*  •  0^ 


1  -i[(Ysint)2-2(l-X-jZ)2] 

(12) 


Y2sintcost[ao(a2a5-b2b5)  -  V^bj+b^)] j 


I 


where  aQ,  bQ  and  all  following  aj,  bj  are  defined  in  the  List  of  Useful 
Formulae.  The  partial  derivatives  of  the  real  part  of  the  phase  re¬ 
fractive  index  with  respect  to  the  spatial  coordinates  (i  =  r,  0,  or  9  ) 
are  similarly  obtained  by  use  of  the  relationship 


Bn  _8n_  _8X  +  .8n_  8_Y+  8_n_8Z  +  Bn  8t 
8i  s  8X  8i  8  Y  8  i  8Z  8i  TfiF  8i 


(13) 
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where 


•'•  ■  |»„  ['2X-‘I  -  <*4*4  -  b4b5»]  *hlZ*  *4b4  *  ^l}  ' 


5  4  5  4 


It  .  RelM  ,  Re '  - 


u  ~ 

3  Y  3  Y 


2 

Re  i  Y  [  (sin^)2  -  i[Y2sin4t+  2(l-X-jZ)2cos2^  ]  ] 

\  MD  L  o  -1 


<  an 
"  a  y 


|  MD  U'  S  - 

Y  |(aoa5  ■  W  [  a6  ’  (8in  )2  ]  "  Vaob5  +  boa5^ 


(15) 


3n 

az 


cos 


an 

Tz 


=  [bo(X-a5a?  +  b5b?)  -  ao(b5a?  +  a^)] 


*)]| 

(16) 


and  where  Re  and  Im  represent,  respectively,  the  real  and  imaginary 
part  of  the  complex  expression.  The  partial  derivative  of  the  angle 
^  (which  is  the  angle  defined  by  the  inner  product  of  the  normalized 
geomagnetic  field  vector  ^  and  the  wave  normal  vector  t  ),  with 
respect  to  the  spatial  coordinates  r,  9,  and  9  ,  measures  the 
change  in  spatial  direction  of  the  earth's  magnetic  field  since  the 
calculation  is  made  holding  8  constant.  The  partial  derivatives 
8X/8i,  3Y/ai  and  8Z/3i  are  obtainable  from  the  analytical  ex¬ 
pressions  for  the  spatial  variation  of  the  electron  density,  geomag¬ 
netic  field  and  collision  frequency.  Examples  of  these  will  be  con¬ 
sidered  later. 
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In  addition  to  Equations  1  to  6,  which  define  the  spatial  ray  path, 
it  is  usually  desirable  to  calculate  the  optical  path  length  s,  the 
time  of  travel  T,  as  well  as,  the  one  way  absorption  A,  of  the 
energy  of  an  electromagnetic  pulse.  The  equation  describing  the 
differential  optical  path  is  given  by 


In  determining  the  time  of  travel  T,  a  distinction  must  be  made 
between  two  electromagnetic  wave  velocities.  The  phase  velocity, 
Vp  =  c/p,  is  defined  as  the  spatial  velocity  with  which  a  point  of 
constant  phase  moves.  Group  velocity,  Vg  =  doo/d  (a)/vp),  is  the 
spatial  velocity  of  electromagnetic  energy  travel;  put  into  other 
words,  it  is  the  velocity  of  a  "Maxwell  Demon"  who  remains  at 
the  same  point  on  the  envelope  of  the  advancing  wave.  From  these 
two  definitions  it  can  be  easily  shown  that  the  time  (in  seconds)  of 
energy  pulse  travel  can  be  written  as: 


dT 

dT 


+ 


(0  9P  "I 
p  9a>  J  ’ 


(18) 


where 

9p  JIM 

9(0  =  Re  9(0  = 

X(2X  +  jZ)  1-  (M2-l)  [2-2jZ-jZX  +  Z{Y^os^>  (l-X-jZ)(l+X) 

S&  *  *  5  [»„<2x2  ‘  *5*8  *  »5V  *  bo<XZ  +  b5*8  *  b8*5»  ]  <19> 


-  Re 


1 


MDCO 
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The  one-way  absorption,  A  (in  nepers),  suffered  by  the  energy  of 
an  electromagnetic  pulse  is  determined  by 


dA 

dx 


-  —  —  A 
c  H 


(20) 


where  k  (which  is  proportional  to  the  spatial  conductivity)  is  the 
absorption  of  the  wave  per  unit  length  of  path. 

The  solution  of  this  set  of  first  order  partial  differential  equations 
will  describe  the  propagation  characteristics  of  an  electromagnetic 
wave  in  a  heterogeneous  anisotropic  medium. 
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B.  LIST  OF  USEFUL  FORMULAE  FOR  RAY  TRACING 


As  an  aid  for  the  computer  solution  of  these  differential  equations 
the  following  list  of  formulae  are  found  to  be  very  useful.  As 
before.  Re  and  Im  ’-espectively  represent  the  real  and  imaginary 
parts  of  the  complex  quantity. 

ReS  =  S.  -  R_  cos  0  (21) 

1  b  5 


ImS  =  S,  =  R  sin  0 

^  b  O 


(22) 


Rs  =  I  [  (Y  sint)4  +  (2Y  cost)2  [  (1  -  X)2  -  Z2]] 


V  2 

+  [(2Y  cost)2  [  2(1  -  X)Z]  J 


1/4 


«  1  -1 
9S  =  2  tan 


(2Y  cost)2  [  2(X  -  1)Z  ] 


(Y  sint)4  +  <2Y  cost)2  [(1  -  X)2  -  Z2] 


ReD 


ImD 


=  dj  =  |z  [(1  -  X)-  Z2  ]  -  (Y  sint)2  +  Sjj 

=  d2  =  [S2  '  2Z<2-X>] 


(23) 


(24) 


(25a) 

(25b) 


ReM  =  m  =  M  =  R  cos  0  (26) 

1  mm 

ImM  =  m  =  -  K  =  R  sin  0  (27) 

2mm 
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Rm  {  1  - 


2X  [(1  -  X)dj  -  Zd2]\‘ 


2  2 
dl  +d2 


2X [ Zd 2  +  (1  -  X)dJ\21 


2  2 
dl  +d2 


2J 


1/4 

}  (28) 


0  =  J tan"1  ( 

m  i.  1 


2X‘[zd1  +  (1  -  X)d£  ] 
d,2  +d22  -  2X[(1  -  X)d,  -  ZdJ 


(29) 


«mldl  -  m2d2) 


<mldl  -  m2d2)2  +  (mld2  +  m2dl)2 


(30) 


b  = 
o 


<mld2  +  m2dl) 


(m1dl  -  m2d2)2  +  (m1d2  +  m^)2 


(31) 


ai  = 


2  [(1  -  X)2  -  Z2]  -  (Y  sin*)2 


(32) 


bj  =  4Z(  1  -  X) 


(33) 


,  .  (aiSl  "  W  1 

1  2  2  J 

S1  +S2 


(34) 
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_  2  ^  _  2 

S1  +  S2 


(35) 


a4  =  < 


b4  = 


1  ♦  2-<Y2c°e’t'|22  [s,(i  -X)  -  ZS2  ] 


s  +  s 
1  2 


Z  +  -2I^°^)22  [s2(1  -  X)  +  zsji 


S1  +S2 


2X  [{1  -  X)dj  -  ZdJ 


2  2 
dl  +d2 


=  (1  +  m22  -  m j2) 


2xfzdj  +  (1  -  X)dJ 


b5  =  d  2  2 

dl  +  d2 


=  2m. 


(36) 


(37) 


(38) 


(39) 


St  |(Y8in2|)2  +  2cos2t[(l  -  X)2  -  Z2]j  -  S2  [(2cob  ty2Z(l  -  X)] 


2  2 
S  +  S 
1  2 


(40) 


b6  = 


S2  |(Y»in2^  )2  +  2cos2i|r[(l  -  X)2  -  Z2  +  Sj  [(2cos*)ZZ(l  -  X)] 


S  2  +  s  2 

S1  5  2 


(41) 
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a?  =  (  (2  -  X) 


+  ^COSt  [s  (1  -  X)  -  s  z] 
Sj  +  s  *  L  1  2  Jl 


(42) 


b?=  2Z+1ML  [siZ  +  S2(1-X)] 


S1  +S2 


(43) 


21  +TFT 71  [n-x2)Sl-s2z(i  +  X)] 


S1  +S2 


(44) 


b8  = 


Z(2  +  X)  +  [SlZ(l  +  X)  +  S2(l  -  x2)  ]  l  (45) 


S1  S2 


I 

I 


’4 
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Nomenclature  Used  in  Ray  Trace  Equations 

r,  0,  cp  spatial  coordinates  of  a  spherical  system 

<?  wave  normal  or  refractive  index  vector 

CT  ,  OQ.CTq)  vector  components  (H) 

X  time  of  phase  travel  (in  units  of  length) 

JU  real  part  of  complex  phase  refractive  index 

f  electromagnetic  wave  frequency 

c  velocity  of  light  in  vacuum 

mass  and  charge  of  an  electron 
real  part  of  the  complex  expression 
imaginary  part  of  the  complex  expression 
complex  phase  refractive  index 


m,  e 

Re 

Im 

M 

X 


(0 

F 

TO 

N 

Y 

TO 


V 

* 


scalar  quantity  representing  normalized  electron 
density 

plasma  frequency  at  the  spatial  point 

angular  wave  frequency 

electron  density  at  spatial  point 

normalized  magnitude  of  the  earth's  magnetic 
field  vector  ? 

magnitude  of  the  gyromagnetic  frequency  of  the 
electron  about  the  earth's  magnetic  field 

normalized  collision  frequency  (v/to) 

collision  frequency  at  spatial  point  in  collisions/ sec. 

angle  defined  by  inner  product  of  vectors  ~0  and  ? 
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v  imaginary  part  of  the  complex  phase  refractive 

index 

k  absorption  coefficient  of  the  wave  per  unit  length 

of  path  (proportional  to  the  conductivity  of  the 
medium) 

i  represents  each  of  the  spatial  spherical  coordinates 

r,  9,9 

T  time  of  travel  (seconds) 

s  optical  path  length  (km) 

A  one-way  absorption  (nepers) 

v  spatial  velocity  with  which  a  point  of  constant 

^  phase  moves 

v  group  velocity  -  spatial  velocity  of  electromag- 

®  netic  energy  travel 
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C.  COORDINATE  TRANSFORMATION 

If  one  takes  three  orthogonal  planes  intersecting  at  a  point,  one 
knows  that  the  position  of  any  point  S  in  space  is  uniquely  determined 
by  the  three  perpendiculars  from  S  on  these  planes,  each  with  its 
proper  sign.  However,  the  problem  of  selecting  the  most  useful 
orientation  of  such  an  orthogonal  system  is  difficult  since  the  use¬ 
fulness  of  a  coordinate  system  partially  depends  on  the  problem 
definition  and  the  application  of  its  solution.  This  ray  trace  program 
is  designed  as  a  sub- set  of  a  much  larger  computer  programming 
effort?  where  the  earth's  geomagnetic  field  plays  an  important  part. 

To  minimize  the  number  of  computer  transformations  in  the  design 
of  the  over-all  program,  an  earth  centered  spherical  coordinate 
system  (r,  9,  <p)  was  chosen,  whose  z-component  is  coincident  with 
the  magnetic  dipole  axis. 

This  selection  permits  the  application  of  the  computer  program  to  a 
great  many  studies  of  ray  path  problems  because  it  accounts  for  the 
earth  curvature  and  accepts  for  solution  any  electromagnetic  radiating 
source  whose  transmitter  location  specifications  of  elevation,  E, 
and  azimuth.  A,  angles,  as  well  as,  geographic  latitude,  9  geo¬ 
graphic  longitude,  A^,  and  position  with  respect  to  the  surface  of 
the  earth  are  known.  Because  the  usefulness  of  this  computer  program 
can  be  extended  by  modification  to  other  coordinate  systems,  as  for 
example,  an  earth  centered  geographic  system,  or  a  radar  coordinate 
system,  the  necessary  coordinate  transformations  from  the  radar  to 
the  earth  centered  geomagnetic  coordinate  system  will  be  described 
in  detail. 

For  the  discussion  of  this  coordinate  transformation,  it  is  assumed 
that  the  electromagnetic  wave  transmitter  is  earth  bound  (i.  e. ,  fixed 
to  the  surface  of  the  earth)  at  a  geographic  latitude  and  a  geographic 
longitude  Ar.  It  is  assumed  that  the  radar  is  positioned  so  that  the 
transmitting  direction  is  described  by  the  elevation  angle,  E,  with 
respect  to  the  tangent  plane  to  the  earth  surface  at  the  radar  location, 
and  an  azimuth  angle.  A,  measured  from  the  radar  coordinate  that 
is  tangent  to  a  great  circle  passing  through  the  north  geographic 
pole.  The  azimuth  angle  is  plus  when  measured  counter-clockwise 
from  the  coordinate  axis,  £,  whose  positive  direction  points  in  the 
direction  of  geographic  north.  It  is  further  assumed  that  S  is  a 
spatial  point  on  the  non-deviated  portion  of  the  ray  a  distance,  R,  from 
the  electromagnetic  wave  transmitter.  This  is  the  spatial  starting 
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point  at  which  the  numerical  methods  necessary  for  solution  of  the 
differential  equations,  must  be  initialized. 

A  transformation  is  required  from  the  spherical  radar  coordinate 
system  to  the  magnetic  coordinate  system  whose  origin  is  at  the 
center  of  the  earth. 

The  necessary  matrices  which  are  required  to  transform  R,  A,  E 
coordinates  to  r,  9,  9  coordinates  can  be  arrived  at  by  a  series  of 
simple  matrix  transformations. 

1.  Let  E,  x]  ,  C  be  a  set  of  orthogonal  coordinates  with  origin 
on  the  surface  of  the  earth.  Let  e-axis  be  perpendicular  to  the 
earth's  surface  while  C  is  directed  (geographically)  northward 
and  q  to  the  east.  As  stated  above  R  is  the  slant  range;  E  is  the 
elevation  angle;  and  A  is  the  azimuth  angle.  Hence  going  from 
R,  E,  A  to  e,  q  ,  C 
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2.  Let  xj,  yj,  z\  equal  an  orthogonal  coordinate  system  with 
origin  on  the  earth's  axis  of  rotation.  The  Xj-axis  is  in  the 
latitude  plane  of  the  radar  site  and  passes  through  the  radar 
site.  The  zj-axis  is  coincident  with  the  north  geographic  co¬ 
ordinate  w.  A  translation  and  rotation  is  required  in  going 
from  e.  q.  t  to  xj,  yj,  Zj. 


where  rQ  equals  the  earth's  radius. 


I 
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3.  Let  X£ ,  Y2*  ^2  etlual  the  orthogonal  coordinate  system  with 
origin  at  the  earth's  center.  Let  the  X2*&xis  be  parallel  to  the 
x^-axis  and  22  coincide  with  w,  hence  also  with  Zj.  Then  going 
from  x^,  yj,  zj  to  X2,  Y2»  z2  by  translation 


(48) 


4.  Let  Xj.  yy  Zj  represent  an  orthogonal  coordinate  system 

with  origin  at  the  earth's  center  such  that  the  Xj-axis  intersects 
the  zero  degree  longitudinal  geomagnetic  meridian  while  the 
z 3-axis  coincides  with  w.  Hence  going  from  X2,  yy  z2  to 
x3*  Y3»  z3  by  rotation  about  the  Z2~axis  yields. 


where  as  before  and  represent  the  geographic  longitude 
and  latitude  of  the  geomagnetic  north  pole  M. 
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Figure  1e.  Geometry  of  Coordinate  Transformation 


5.  Let  x,  y,  z  represent  an  orthogonal  coordinate  system  with 
origin  at  the  earth's  center.  Let  the  x-axis  pass  through  the 
great  circle  connecting  the  geographic  and  geomagnetic  poles 
while  the  z-axis  passes  through  the  geomagnetic  pole  M,  In 
going  from  X3,  y^,  Zj  to  x,  y,  z  by  a  rotation  one  obtains 


M 

1  .in*M 

0 

"co*^M 

y 

s 

0 

1 

0 

n 

i  C°’fM 

0 

X  \ 

1  X 

1 

3 

y 

■  <v 

y3 

z 

1  z 

1 

\  3 

(50) 
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Figure  Id.  Geometry  of  Coordinate  Transformation 


Hence  by  matrix  multiplication  one  can  transform  from  radar  co¬ 
ordinates  R,  E,  A  to  the  earth  centered  coordinates  x,  y,  *  where 
the  magnetic  dipole  axis  of  the  earth  coincides  with  the  z-axis.  From 
this  coordinate  system  one  can  simply  transform  to  the  desired 
spherical  coordinate  system  r,  9,  tp- 
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'ij 


<V 


gij  = 


ij)  <v  (bij)  ■  (fij) (b 

u1 

(52) 

sinq>Mcos(AM-AR) 

_cos9m\ 

-8in(AM-  V 

co»(Am-Ar) 

0 

(53) 

co8<Pmcos(Am-Ar) 

co^M.m(AM  -Aft) 

Sin,PM  j 

/ 8 1 1  g12 

gl3\ 

g21  g22 

g23 

(54) 

^831  g32 

833 

gu  =  (co8fR  8in<PM  C08(AM-  A^  -  C08<PM  8in<PR] 

8lZ  =  8infM  8in(AM  •  V 

813  =  l'Sin<PR  8in<PM  C08(  AM"AR)  ’  COS<PR  C°B9m) 
821  =  -8in(AM  -AR)c08<f>R 

«22  =  C08(Am'V 

«23  =  sin^R  8in^  AM  "  A  R> 

831  =  |COs9R  C08?MC08(AM  '  V  +  8in9M  8in9R) 

832  =  C°8*M8in(AM‘AR> 

833  =  (-Sin?R  C089MC08(AM-AR)  +  “in*M  CO'*r) 
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' 

sin  E  ^ 

sin  9  cos<A 

y 

■  RV 

cos  E  sin  A 

+  r»  <•«>  ■  r 

sin  9  sin9 

\z 

\  cos  E  cos  A 

\  cos  0  / 

By  use  of  these  matrix  transformations  the  Cartesian  coordinates, 

(x,  y,  z),  and  from  them  the  spherical  coordinates,  (r,  0,  9),  of  the 
earth  centered  geomagnetic  coordinate  system  can  be  determined  for 
the  spatial  starting  point  S  and  the  earth  bound  transmitter  R.  These 
can  be  expressed  as: 


x 


S 


+ 


[cos(AM-AR)  sin9M  cos<pR  -  cos<PM  sinq>R  ]  (R  sin  E  +  rj 
jsin(Aj^  -  Aj^)  sintp^j  R  cos  E  sin  A  (56) 

^cos(AM  -  Ar)  sin«PM  sin<PR  +  cos«p^  cos<PRJ  R  cos  E  cos  A 


yS  *  [’  8in(AM  '  V  COB*r\  (R  E  +  ro)  +  [C08(AM-  Vl  R  cos  E  sin  A 
+  ^sin(A^  -AR)  s in  9 R|  R  cos  E  cos  A  (57) 

*S  *  [CO,(AM  •  V  CO'TMCO,lfR  *  ,in»M,in,J 

+  |sin(AM  *  AR)cos«PM]  R  COS  E  sin  a 

-  h"AM  ‘  VC<>,,M  ,in,,R  ‘  ,in9M  co,,’r] 


(R  sin  E  +  r  ) 
o 

(58) 

R  cos  E  cos  A 
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When  R  =  0,  that  is,  for  a  point  on  the  surface  of  the  earth, 


XR  *  ro  [COs(AM  -  AR>  8in*Mc°8<PR  -  cosq)M8in<pR]  (59) 

yR  =  -  rQ  sin(AM  -AR)  cos<pR  (60) 

zR  =  rQ  [cos(AM  -  AR)cos«PMcos<pR  +  sin^sin^]  (61) 


From  simple  trigonometric  considerations  (Figure  2a)  it  can  be 
shown  that  the  radar  slant  range,  R,  measured  from  the  transmitter 
to  the  spatial  starting  point  S  is  given  by 


R  = 


r  sin  E  + 
o 


V 


2 

r  cob  E 
o 


(62) 


where  1^  is  the  vertical  height  of  the  starting  point  above  its  pro¬ 
jection,  (point  P)  on  the  surface  of  the  earth. 

Equations  1  through  6  point  out  that  in  addition  to  these  transforma¬ 
tions,  the  components  of  the  directed  normal  to  the  phase  fronts, 

3  ,  at  the  starting  point  S  are  to  be  determined  in  this  coordinate 
system.  From  spherical  and  plane  trigonometric  considerations, 
(Figure  2),  it  can  be  shown  that  these  components  are  given  by 


a  cos  e 

(63) 

0  sin  e  cos  a 

(64) 

-  0  sin  e  sin  a 

(65) 
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Angle  e  can  be  evaluated  directly  by  employing  the  law  of  sines. 
This  yields 


e 


.  -1 

sin 


'  r  cos  E 
o _ 

r  +  h_ 
i  o  S 


(66) 


Angle  a  is  the  geomagnetic  bearing  angle  (Figure  2)  measured  posi¬ 
tive  in  a  clockwise  direction  from  geomagnetic  north.  By  use  of 
spherical  trigonometry  it  is  expressible  by 


a 


tan 


-1 


8in(*S  ~  $R 
cos  0gSin9^cos($g 


)sin0R 

-  $R)-  sin0scos0R 


(67) 


where  $R,  0R,  and  $g,  0g  are  the  geomagnetic  longitudes  and  co¬ 
latitudes,  respectively,  of  the  radar  transmitter,  R,  and  the  spatial 
starting  point,  S.  The  geomagnetic  angles  are  obtained  from  the 
following  expressions 


-1 

-1 

^g2l\ 

=  tan 

(X*J 

=  tan 

l®11/ 

=  tan 


-1 


"in,AM  *  AR,CO*,'r 


LCO,,AM-V,in  V^R  *  '°<Tm»‘»*r 


(68) 


0R  =  co. 


=  COS 


■1  I  R  \  -i  ,  . 

=  cos  (g,,) 
lrR  |  '*31' 


1 


CO,<  V  AR)CO”  MCO,<fR  *  si"<PM.in0  R 


(69) 
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t 


S 


(70) 


0 


S 


-1 

cos 


z 


S 


r 

o 


(71) 


One  additional  useful  expression  can  be  obtained  from  these  algebraic 
relations.  The  parameter  is  the  angle  ^  at  the  spatial  starting 
point  S  which  is  defined  by  the  inner  product  of  the  magnetic  field 
vector  and  the  wave  normal.  By  the  application  of  the  sine  and  cosine 
laws  to  »he  geometry  of  Figure  2-a,  it  can  be  shown  that 


*  = 


cos 


1  [-(cos 


e  sin  I  +  sin  e  cos  I  cos  « 


>] 


(72) 


where  angle  I  is  the  magnetic  inclination  angle.  The  inclination 
angle^  is  only  a  function  of  the  geomagnetic  latitude  at  the  parti¬ 
cular  point  in  question. 


I  =  tan 


'[ 


2  cotan  8, 


(73) 


Expressions  arising  from  the  inverse  coordinate  transformation, 
that  is,  transformation  from  the  geomagnetic  coordinates  (r,  0,«p) 
to  the  radar  coordinates  R,  E,  A  can  be  easily  developed  from 
these  formulae.  Although  used  in  the  computer  program  they  will 
not  be  presented  here. 
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Figure  2.  Geometry  for  Storting  Point  and  Geomagnetic  Coordinate  System 
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Nomenclature  Used  in  Coordinate  Transformation 


r,  0,<P 


E 

A 

<p 

R,  M 

ar,  M 
S 

R 


e»  n-C 


spatial  point  in  an  earth  centered  spherical  co¬ 
ordinate  system 

radar  elevation  angle 

radar  azimuth  angle 

geographic  latitude  of  point  R  or  M,  respectively 

geographic  longitude  of  point  R  or  M,  respectively 

spatial  starting  point  on  nondeviated  portion  of  ray 

distance  from  electromagnetic  wave  transmitter 
to  starting  point 

set  of  orthogonal  coordinates  with  origin  on  the 
surface  of  the  earth  at  radar  site 


Xl'  yr  Z1  orthogonal  coordinate  system  with  origin  on  the 

earth's  axis  of  rotation 


x2'  y2‘  z2  orthogonal  coordinate  system  with  origin  at 

earth's  center 


3’ 


r3*Z3 


orthogonal  coordinate  system  with  origin  at  earth's 
center 


w  z  component  of  the  geographic  coordinate  system 

(u,  v,  w) 

hg  height  of  starting  point  above  its  projection  on 

the  surface  of  the  earth 


r  radius  of  earth 

o 


a 


geomagnetic  bearing  angle 


^R’ 


geomagnetic  longitudes  of  points  R  and  S 


29 


0R'eS 


°r*  V°«P 


geomagnetic  co-latitudes  of  points  R  and  S 

physical  components  of  a  vector  of  length  P,  that 
is  directed  normal  to  the  phase  front 


* 


angle  between  magnetic  field  vector  and  the  wave 
normal 


I 


angle  of  magnetic  inclination 
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D.  MODEL  IONOSPHERE 

As  shown  under  Computational  Procedure,  the  refractive  index  M 
and  its  spatial  derivatives  are  dependent  on  the  normalized  density, 
X,  and  its  spatial  gradients,  8X/9r,  9X/9  9,  and  9X/9?  .  For 
an  evaluation  of  these  quantities  an  analytic  model  ionosphere  can 
be  chosen.  One  such  ionospheric  model  that  was  found  useful,  is  a 
spherical  electron  distribution  as  measured  from  a  point  in  space. 

A  reason  for  its  selection  is  presented  under  Computational  Results. 
Other  uses  as  well  as  other  ionospheric  models  are  covered  else¬ 
where^*  1 1 . 


Let  r^,  0^,  Cpfc  represent  the  spatial  location,  B,  of  the  center  of 
the  selected  spherical  ionosphere  in  the  geomagnetic  spherical  co¬ 
ordinate  system  (r,  0,  ?  ).  The  values  of  these  coordinates  are 
obtainable  from  the  specified  geographic  latitude,  longitude  and 
height  above  the  earth  surface  of  point  B,  in  an  analogous  procedure 
as  described  for  the  transformation  of  coordinates  of  the  spatial 
starting  point  S.  Then  the  electron  density  at  a  spatial  point  r,  0,? 
for  an  assumed  spherical  ionosphere  can  be  written  as: 

N«r-  »■  *•  v  v  v  ,74) 


where 


((rb  sin  0fa  cos?b  -  r  sin  0  cos?)  +  (rbsin0bsin?b  -  r  sinOsin?) 


1/2 


1/2 


+  (r.  cos  0,  -  r  cos  0) 

D  D 


=  <XP2  +  YP2  +  ZP2 


(75) 


For  this  discussion  A  and  n  are  appropriately  chosen  constants  which 
give  the  desired  electron  density  N  (electrons/cc).  By  use  of 
Equations  74  and  75,  it  is  easily  shown  that  the  spatial  electron 
density  gradients  can  be  expressed  in  the  following  manner. 
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3N 

9r 


nA 


„n  +  2 


XP  sin  0  cos  9  +  YP  sin  8  sinf  +  ZP  cos  0 


(76) 


8N 

5F 


nAr 
n  +  2 


t 


XP  cos  0  cos<p  +  YP  cos  0  sin  9  -  ZP  sin  0 


8N 

89 


nAr 
-n  +  2 


-  XP  sin  0  sin9  +  YP  sin  0  cos9 


(77) 


(78) 


Some  results  obtained  by  use  of  such  a  spherical  ionospheric  model 
will  be  discussed  later. 
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E.  MODEL  OF  EARTH'S  MAGNETIC  FIELD 


Because  magneto-ionic  effects  on  the  propagation  of  electromagnetic 
waves  through  an  ionized  medium  are  taken  into  account  in  the 
derivation  of  the  equations  under  Computational  Procedure,  it  is 
necessary  to  specify  the  normalized  external  magnetic  field  of  the 
earth,  Y,  its  components  Y  ,  Yq,  Y«p  and  its  spatial  derivatives 
3Y/9r,  0Y/00,  0Y/09.  It  is  known  that  the  earth's  magnetic 
field  can  be  approximated  b/  an  earth  centered  magnetic  dipole  with 
its  axis  displaced  such  that  the  geographic  longitude  =  70.  1°W 
and  the  geographic  latitude  =  78,6°N.  The  magnetic  potential, 

V,  at  a  distant  point  from  such  a  dipole  is  related  to  the  magnetic 
moment, ,  by  the  expression 


V(r,e)  =  .T&Sfl 
r 


(Y  r  ■*)  cos  8 
e  o 

_ 

r 


(79) 


where  r,  0  are  the  geomagnetic  coordinates  of  the  spatial  point  ir¬ 
respective  of  the  coordinate  9  and  as  before,  rQ  =  radius  of  the 
earth.  In  this  equation  Ye  is  the  magnitude  of  the  normalized 
magnetic  field  at  the  earth's  surface  on  the  magnetic  equator.  By 
use  of  this  algebraic  equation  all  the  desired  quantities  can  be 
derived.  They  are 


Y 


Y  (— )3  (1  +  3  cos2  0)1/2 
e  r 


(80) 


where,  as  previously  defined,  Y  is  the  normalized  magnitude  of  the 
earth's  magnetic  field  vector  ?  =  (eH/mcto)  a  CDc/co 


Y 

r 


2  Y  (— ^  cos  0 
e  r 


1  * 

—  tan  9 
4 


Y 


0 


Y  (— )3  sin  0 
e  r 


2 


Y  tan  0 
r 


(81) 


(82) 
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3Y 

=  5c r  =  0 


(83) 


3Y_  3Y; 

3r  r 


(84) 


3Y 

ae 


3Y  sin  0  cos  0 


1  +  3  cos  0 


(85) 
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F .  MODEL  OF  ATMOSPHERIC  COLLISION  FREQUENCY 


For  some  of  the  trial  calculations  the  atmospheric  collision  fre¬ 
quency  was  found  from  assumed  exponential  variations  of  collision 
frequency  with  height.  The  atmosphere  was  radially  stratified 
and  an  approximate  exponential  equation  was  curve-fitted  to  mea¬ 
sured  experimental  data  for  each  stratified  rejron.  Hence,  for 
each  region  the  following  relations  were  used  t  obtain  Z  and 
9  Z /  9r,  9Z/  90,  9Z/8cp  that  are  required  by  the  ray  trace  equations. 


Z 


CO 


-b(r  -  r  ) 
ae  o 


(86) 


az 

3r 


-bZ 


(87) 


az  az 
ee  =  a«p 


o 


(88) 


The  dependence  of  collision  frequency  on  a  localized  temperature 
distribution  and  degree  of  ionization?  complicates  these  simple 
relations.  These  complications  (as  derived  by  D.  Archer)  as  well 
as  their  effects  will  not  be  discussed  at  this  time. 
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G.  COMPUTATIONAL  RESULTS 

The  preceding  equations  are  only  a  summary  of  the  required  set 
which  will  permit  the  detailed  calculation  of  a  ray  path  in  three- 
dimensional  space.  Because  of  this,  it  becomes  clear  that  the  only 
realistic  approach  to  the  solution  of  this  problem  is  the  utilization 
of  computer  techniques,  otherwise,  the  welter  of  data  that  must  be 
handled  through  use  of  numerical  methods,  is  beyond  effective 
human  handling  capacity.  However,  the  development  of  a  computer 
program  which  can  perform  countless  number  of  calculations,  poses 
the  very  difficult  task  of  determining  the  correctness  of  a  computed 
result.  To  simplify  this  "debugging"  task  the  classical  idea  of 
elastic  collision  between  charged  particles  was  borrowed  from  nuclear 
physics.  It  has  been  shown! 0  that  if  the  electron  density  falls  off  as 
the  inverse  square  of  the  distance  from  the  center  of  a  spherical 
electron  distribution  (N  =  various  exact  expressions  can  be 

obtained,  since  the  ray  equations  at  zero  azimuth  are  integrable. 

The  geometry  of  such  a  distribution,  as  well  as,  three  ray  paths 
computed  by  use  of  the  computer  program  are  shown  in  Figure  3.  In 
the  figure  the  center  of  the  sphere  is  located  by  the  fixed  coordinates 
(R0.  0O)  with  respect  to  the  radar.  The  derivations  are  made  in  two 
dimensions,  hence  only  the  two-dimensional  coordinate  system  ( £  ) 
is  used.  The  ray  path  has  an  initial  elevation  angle  Ej  and  its  distance 
of  closest  approach  to  the  center  of  the  refracting  sphere  is  denoted 
by  The  coordinates  of  any  point  on  the  ray  path  with  respect  to 

the  center  of  the  sphere  are  (*£  0)  and  with  respect  to  the  radar  (R,  E). 
Angle  b  represents  the  amount  of  ray  bending  experienced  by  the 
ray  passing  through  the  refractive  medium.  Under  these  conditions 
Archer!®  has  shown  that  the  angle  &  is  given  by 


S  =  u  -  2  [ro  +  n  (^o)  cos-1  (-—-)]  (89) 

o 

The  three  plotted  rays  have  actually  a  small  third  dimensional  com¬ 
ponent.  Because  of  this,  the  accuracy  of  these  plotted  rays  is  approxi¬ 
mately  (+3,  -2)  percent. 

As  shown  by  the  tabulated  results  in  Table  1,  the  computed  deviation 
angle  &c,  arising  from  ray  trace  results,  agrees  very  well  with  the 
calculated  angle,  &,  obtained  by  use  of  Equation  89.  These  computer 
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calculations  were  performed  under  a  large  error  upper  bound  con¬ 
dition  (E  =  10* 3,  see  description  of  Subroutine  INT).  The  agreement 
is  improved  by  a  variation  of  this  error  condition.  Additional  results 
will  not  be  presented  here.  Presented  elsewhere^  is  the  influence  of 
E  on  computed  results  using  these  numerical  methods  as  applied  to 
the  study  of  ionization-deionization  phenomena. 


Radar  Elevation 

Bending  Angle  -  & 

Bending  Angle  -  &c 

Angle  -  E 

from  Equation  89 

Ray  Trace  Program 

Degrees 

Degrees 

Degrees 

60 

19.  5 

20 

70 

35.  2 

37 

80 

77.  5 

78 

Table  1 

Comparison  of  Total  Ray  Bending  Angle 


As  an  additional  illustration,  refractive  errors  through  a  particular 
spherical  ionized  region  have  been  computed  in  detail  to  illustrate 
the  concepts  discussed  and  to  indicate  the  kinds  of  refractive  errors 
that  could  arise.  The  electron  density  contours  of  this  ionization 
model  are  defined  by  Equation  74  where  A  =  10^3  and  n  =  12.  The 
distance  from  the  center  of  the  spherical  ionosphere  is  measured  in 
kilometers.  The  center  of  the  model  is  located  at  an  elevation  angle, 
E,  of  30  degrees,  zero  degree  azimuth  angle,  and  564  kilometers 
slant  range  as  measured  from  the  radar  site. 

Figure  4  shows  the  relation  between  the  radar  and  ionization  model 
in  the  plane  of  zero  azimuth.  Also  shown  are  the  ray  paths  for  rays 
leaving  the  radar  at  several  elevation  angles.  A  frequency  of  one 
kilomegacycle  was  used  in  determining  the  refraction  of  the  electro¬ 
magnetic  wave  propagation  vector.  Because  the  electron  density 
increases  rapidly  near  the  center  of  the  model,  there  is  significant 
bending  of  the  ray  path. 
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The  refraction  becomes  so  severe  as  the  elevation  angle  of  radar 
rays  approach  the  elevation  angle  of  the  ionization  center  with  respect 
to  the  radar,  that  there  is  a  region  (shown  by  half  tones)  into  which 
no  radar  ray  penetrates,  hence,  radar  "blackout"  is  achieved.  In 
three  dimensions  this  blackout  region  is  a  cone  in  which  a  target  is 
shielded  from  the  radar.  Because  rays  near  this  region  intersect 
each  other,  two  elevation  angle  paths  to  the  same  target  exist,  so 
multiple  targets  may  be  visible. 

If  the  ray  path  is  not  in  the  zero  azimuth  plane,  the  amount  of  elevation 
error,  or  azimuth  error,  is  a  function  of  the  location  of  the  target. 

The  elevation  and  azimuth  errors  for  a  target  located  at  a  slant 
range  of  1200  kilometers  have  been  computed  as  a  function  of  radar 
elevation  and  azimuth  angles.  These  are  summarized  in  Figure  5, 
in  which  contours  of  constant  elevation  error,  Ae,  in  one  quadrant 
and  constant  azimuth  error,  A  A,  in  another  quadrant  as  a  function 
of  the  ray  direction  at  tne  radar  site  are  given.  The  contours  have 
been  terminated  at  a  total  bearing  error  of  about  10  degrees.  Due  to 
symmetry  the  errors  in  the  other  quadrants  are  just  the  mirror  image 
of  the  quadrants  shown. 
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R  distance  between  the  radar  and  center  of  the  sphere 

o 

Yq  apparent  bearing  angle 
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SECTION  III 

COMPUTER  PROGRAM  FOR  THREE-DIMENSIONAL  RAY-TRACING 


Figure  6  schematically  describes  the  computer  program  that  was 
developed  for  three-dimensional  ray  tracing.  As  illustrated,  the 
computer  program  is  a  composite  of  a  group  of  subprograms. 

Because  each  subroutine  is  an  entity  in  itself,  the  improvement  of 
the  entire  program  can  be  performed  by  the  variation  of  each  sub¬ 
program. 

For  the  creation  of  this  program  the  FORTRAN  language^  was 
used  wherever  possible.  FORTRAN  is  an  automatic  coding. system 
for  the  IBM-704/709/7090  Data  Processing  Computer  System  that 
was  designed  for  scientific  application.  Although  there  are  limita¬ 
tions  to  FORTRAN,  nevertheless,  1)  it  is  at  present  the  only  language 
for  scientific  use,  that  is  accepted  by  most  existing  large  computer 
systems,  and  2)  it  is  simple  and  therefore  without  much  effort, 
permits  the  elimination  of  the  programmer,  thus  leaving  the  design 
of  logical  computer  decisions,  to  the  formulator  of  the  scientific 
problem.  The  program  has  been  written  to  operate  "in  or  out"  of 
the  FORTRAN  MONITOR  CONTROL  SYSTEM. 

Except  for  the  RINDEX  subroutine  the  program  has  been  divided 
into  small,  simple  Functions  and  Subroutines  to  facilitate  under¬ 
standing.  In  the  development  of  the  program,  concentration  was 
mainly  on  obtaining  a  correct  working  program,  as  soon  as  possible, 
and  not  on  optimization  or  clarity  of  output  results.  These  tasks 
are  left  for  future  development. 

The  computer  program  consists  of  the  following  parts: 


1) 

Main  Program 

RAY  TRACE 

2) 

Function 

SLANTR 

3) 

Function 

QATAN 

4) 

Function 

ARCOS 

41 


RM  61TMP-32 


5) 

Subroutine 

COORD 

6) 

Subroutine 

DAUX 

7) 

Subroutines 

INT  and  INTM 

8) 

Subroutine 

RINDEX 

9) 

Subroutine 

ELECTX 

10) 

Subroutine 

BIGR 

11) 

Subroutine 

MAGY 

12) 

Subroutine 

COLFRZ 

13) 

Subroutine 

RCOORD 

14) 

Subroutine 

OUTONE 

15) 

Subroutine 

OUTPUT 

These  functions  and  subroutines  are  used  to  obtain  numerical  values 
for  those  variables  which  cannot  be  defined  by  only  one  arithmetic 
statement.  In  addition  to  these  subprograms  certain  statements  in 
the  FORTRAN  language  cause  the  inclusion  in  the  object  program  of 
the  necessary  input  and  output  routines,  as  well  as,  various  library 
functions  and  subroutines  in  relocatable  binary  form  that  are  avail¬ 
able  on  the  FORTRAN  MASTER  LIBRARY  TAPE.  The  names  and 
locations  of  these  necessary  routines  are  given  in  the  "storage  map" 
of  the  arrangement  of  storage  location  in  the  object  program  that  is 
compiled  from  a  FORTRAN  source  program.  These  "maps"  follow 
the  listings  of  each  source  program.  These  added  routines  will  not 
be  discussed  in  this  report.  Following  a  brief  description  of  the 
function  of  the  main  program  and  its  associated  subprograms,  the 
necessary  input  data  for  a  sample  calculation  is  given  with  a  des¬ 
cription  of  the  output.  Some  of  the  results  listed  in  this  output  led 
to  the  graphical  results  presented  in  Figures  3,  4  and  5. 
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Figure  4.  Radar  Propagation  Paths  through  Spherically  Ionized  Region 
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Figure  5.  Elevation  and  Azimuth  Error*  for  Propagation  through 
a  Spherical  Model  -  f  *  I  kMc 


44 


DEGREES 


I 

I 


RM  61TMP-32 


OUTPUTS 


nA7“ 

|  TtACINC 
!  Moot  AM 


DIPOLE 


EMPIRICAL 

MAGNETIC 

DATA 


COORDINATE  TRANSFORMATION  -  RAY  SOURCE 
TO  EARTH  CENTERED  GEOMAGNETIC 


NON-DEVIATIVE  RAY  TO 
IONIZED  CLOUO  BOUNDARY 


MAGNETIC  FIELD  AND  ITS  DERIVATIVES 


BEST  KNOWN  ATMOSPHERE 


COLLISION  FREQUENCY 


THREE  DIMENSIONAL  HASELGROVE  EQUATIONS 
FOR  EM  WAVE  PROPAGATION 


DIFFERENTIAL  EQUATION  SOLUTION  BY 
RUNGE-KUTTA  OR  ADAMS-MOULTON 


R INDEX  -  CALCULATION  OF  INDEX  OF 
REFRACTION  AND  ITS  DERIVATIVES 
AND  THE  ABSORPTION  COEFFICIENT 


ElECTX  <ALCULATE  ELECTRON  DENSITY 
•  AND  ITS  DERIVATIVES  AS  CREATED  BY  ALL 
ION  PRODUCTION  MOOES  FROM  DETONATION 


AMBIENT  IONIZATION 


DEIONX  -  RAH  OF  ELECTRON  LOSS 
AS  FUNCTION  OF  SPECIES  PRESENT 


MOTION  -  EXPANSION  OF  POINT  SOURCE 
CLOUD  RISE,  ANO  SHOCK  WAVES 


44a 
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A.  MAIN  PROGRAM  RAY  TRACE 

The  main  program's  function  is  to  act  as  a  master  control  of  the 
logical  flow  necessary  for  the  execution  of  the  numerical  methods. 

In  addition,  it  is  responsible  for  obtaining  the  necessary  data, 
initializing  the  required  starting  conditions,  performing  the  desired 
controlled  printouts  of  computed  results,  and  determining  the  condi¬ 
tion  for  termination  of  the  given  computations.  For  initializing  the 
starting  conditions  the  main  program  requires  the  reading  into 
storage  of  the  following  information  in  the  format  illustrated  under 
INPUT,  Table  4. 


CARD  1  RECORD:  This  can  be  any  desired  infor¬ 

mation  consisting  of  72  alphanumeric 
characters  that  will  serve  to  identify  the 
calculations. 

CARD  2  Contains  the  values  of  ID,  KWIT.  ID  is 

an  integer  that  can  be  used  for  identifying 
the  calculation  if  the  same  CARD  1  is 
used.  KWIT:  On  completion  of  the  calcu¬ 
lations  in  progress,  the  computer  will 
check  if  additional  problems  are  to  be 
performed.  Thus  if  KWIT=59  the  computer 
will  want  to  read  a  new  W  vector 
(CARD  3  —  onward);  if  KWIT=66,  the 
computer  will  want  to  read  a  new  CARD  1, 
CARD  2,  and  new  W  vector  (CARD  3  —  on¬ 
ward);  if  KWIT  equals  any  other  integer 
the  computer  will  PAUSE  44444. 

CARD  3  onward  This  card  and  all  following  cards  des¬ 

cribe  the  value  of  each  component  of  the 
W  vector  that  is  not  zero.  As  shown 
under  INPUT,  the  first  three  columns 
of  the  card  are  for  the  integer  that 
describes  the  W  vector  component.  The 
next  fourteen  columns  of  the  card  are 
for  the  value  of  the  W  vector  component. 
The  number  of  these  cards  is  variable 
since  on  completion  of  one  calculation, 
often  only  one  component  of  the  W  vector 
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CARD  3  (continued) 


LAST  DATA  CARD 


SENSE  SWITCH  1 


SENSE  SWITCH  2 


SENSE  SWITCH  3 


SENSE  SWITCH  4 


SENSE  SWITCH  6 


is  to  be  changed  for  the  next  computa¬ 
tion,  The  W  vector  can  be  read  in  any 
order. 

This  card  follows  the  last  card  descrying 
the  W  vector.  It  is  any  negative  nteger 
listed  in  the  first  three  columns  o  a 
card.  It  transfers  the  computer  o^t  of 
the  read  mode  to  the  location  beginning 
the  ray  trace  calculations. 

The  program  is  designed  to  calculate 
first  the  ordinary  ray  path  and  then  the 
extraordinary  ray  path.  SENSE 
SWITCH  1  DOWN  will  eliminate  the 
calculation  of  the  extraordinary  ray 
path. 

In  DOWN  position  will  permit  the  cal¬ 
culation  of  the  extraordinary  ray  path 
and  eliminate  the  calculation  of  the 
ordinary  ray  path. 

Placing  this  SENSE  SWITCH  3  DOWN 
will  terminate  the  calculation  on  com¬ 
pletion  of  the  ray  path  calculations  in 
progress. 

Placing  SENSE  SWITCH  4  DOWN  will 
cause  the  computer  to  check  if  SENSE 
SWITCH  6  is  DOWN.  If  it  is  down  the 
computer  will  terminate  calculations 
immediately. 

It  is  desirable  to  follow  the  course  of 
any  calculation  on  a  computer.  SENSE 
SWITCH  6  DOWN  will  print  on-line, 
the  total  number  of  numerical  integrations 
completed  up  to  this  point,  integration 
mesh  size,  length  of  independent 
variable  “t ,  height  above  surface  of  the 
earth  (km),  0,  ? (in  degrees),  (Jr,  0g, 

,  P,  * ,  distance  from  the  ion  source 
center  (km),  value  of  the  normalized 
electron  density  X. 
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SENSE  LIGHT  2 


PAUSE  17171 


PAUSE  66666 


PAUSE  44444 


When  SENSE  SWITCHES  1  and  2  are  DOWN, 
then  both  ordinary  and  extraordinary  ray 
paths  are  being  calculated.  When  SENSE 
LIGHT  2  is  ON,  then  the  calculation  is 
determining  the  ordinary  ray  path.  When 
it  is  OFF,  then  the  extraordinary  ray 
path  is  being  calculated. 

If  the  computer  halts  with  this  octal 
number  in  the  address  field  of  the 
STORAGE  REGISTER  it  signifies  that 
SENSE  SWITCHES  1  and  2  are  in  the  UP 
position  and  the  problem  is  undefined. 
SENSE  SWITCHES  1  or  2, or  1  and  2  are 
to  be  placed  DOWN  depending  if  only  the 
ordinary,  the  extraordinary,  or  both  ray 
paths  are  to  be  calculated.  Following  the 
definition  of  the  problem,  pushing  START 
key  will  cause  calculations  to  resume. 

The  computer  halts  with  this  octal  number 
in  the  STORAGE  REGISTER  just  prior 
to  beginning  calculations.  If  the  MONITOR 
system  is  used  it  permits  the  operator 
to  know  when  it  has  left  the  MONITOR 
system  and  the  SENSE  SWITCHES  can  be 
changed  as  needed  by  the  problem. 

The  computer  halts  with  this  octal  number 
in  the  STORAGE  REGISTER  on  completion 
of  all  the  necessary  calculations  specified 
by  the  INPUT  data.  It  permits  the  opera¬ 
tor  to  reset  the  desired  sense  switches 
for  the  MONITOR  system.  Pressing 
START  will  cause  the  computer  to  exit 
from  the  program  to  the  MONITOR  system. 


Table  2  contains  the  nomenclature  that  describes  some  of  the  com¬ 
ponents  of  the  V  vector  and  the  components  of  the  W  vector. 
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V(2) 

independent  variable  X 

V(3) 

initial  step  size  input  Ax 

V  (4) 

radius  from  center  of  earth  r 

V  (5) 

variable  angle  8 

V(6) 

variable  angle  9 

V  (7) 

o 

r 

V(8) 

ae 

V(9) 

V(10) 

optical  path  length  one  way  s 

V{11) 

time  one  way  T 

V(12) 

A  absorption 

V(  13) 

dr/dT 

V(14) 

d0/dx 

V(  1 5) 

d<p/dx 

V(16) 

do  /dx 
r 

V(17) 

dO0/dx 

V(18) 

do?/dx 

V  ( 19) 

ds/dX 

V(20) 

dT/dX 

V(21) 

dA/dX 

Table  2.  Nomenclature  Describing  the  V  and  W  Vectors. 
(Page  1  of  6) 
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W(l) 
W(2) 
W(3) 
W  (4) 
W(  5) 
W(6) 
W(7) 
W  (8) 
W  (9) 
W(10) 

wcii) 

W(12) 

W(13) 

W(  14) 
W(15) 

W(  16) 
W(17) 
W(18) 
W(19) 
W(20) 


refractive  index  H 

imaginary  part  of  complex  phase  refractive  index  X 

radar  transmitter  angular  frequency  to 

3p/8a 

r 

au/ao0 

8u/3r 
3  H/30 
9p/0<p 
3n/9ir 
Ota/Oco 

unassigned  for  this  program 

geographic  longitudinal  angle  A^  of  geomagnetic 
north  pole  measured  east  of  Greenwich  Meridian 
(degrees) 

angle  measured  as  W(13)  in  degrees 

angle  geographic  latitude  of  geomagnetic  north- 
pole  measured  plus  from  geographic  equator  north 

angle  tP  p  geographic  latitude  of  radar  (degrees) 

radar  elevation  angle  E  (degrees) 

radar  azimuth  bearing  angle  angle  A  (degrees) 

rQ  radius  of  the  earth  (km) 

hg  height  of  starting  point  above  surface  of  earth 

Table  2.  Nomenclature  Describing  the  V  and  W  Vectors. 
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RM  61TMP-32 

W(21 ) 

W  (22) 

W  (23) 

W  (24) 

W(25) 

W  (26) 
W(27) 

W  (28) 

W(29) 

W{30) 

W(31) 

W(32) 

W(33) 

W(34) 

W(35) 

W(36) 

W(37) 

W(38) 


angle  of  ionization  source  measured  as  W(15) 
{degrees) 

longitudinal  angle  Ag  of  source  measured  as  W(13) 

hg  height  of  ionization  source  center  above  earth 
surface 

At  initial  mesh  size  of  variable 

Yg  normalized  equator  magnetic  field  on  earth's 
surface  at  the  geomagnetic  equator 

a  constant  determining  collision  frequency 

b  constant  in  exponent  determining  collision 
frequency 

range  (km)  =  distance  from  ionization  source  center 
to  spatial  point  (r,  8,<P  )  = 

cosine  of  angle  makes  with  the  vertical  through 
center  of  the  ionizing  source 

radial  distance  from  earth's  center  to  center  of 
ionizing  source 

x  geomagnetic  coordinate  of  source 

y  geomagnetic  coordinate  of  source 

z  geomagnetic  coordinate  of  source 

A  =  constant  in  A/R°  determining  electron  density 

n  =  exponent  in  A/Rn  equation 

unassigned  for  this  program 

unassigned  for  this  program 

plasma  angular  frequency  cycies/sec 

Nomenclature  Describing  the  V  and  W  Vectors. 
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RM  61 TM  P-32 


W(39) 
W  (40) 
W(41) 


W  (42) 

W(43) 

W  (44) 

W(45) 

W(46) 

W  (47) 

W(48) 

W(49) 

W(50) 


Ng  in  ion  pairs /cc 

maximum  r  =  V(4)  to  be  considered  in  this  calculation 
A1  =  vector  in  INTM  routine 

If  W{41)  =  0  routine  will  use  predictor  corrector  with 
variable  VT(24) 

If  W(41)  =  2  will  use  Runge-Kutta  with  fixed  W(24) 

If  W(41)  =  2  will  use  predictor-corrector  with  fixed 
W  (24) 

If  W(41)  =  1  or  2  then  W(42)  through  W(47)  are  ignored 
but  must  have  some  value. 

If  W(41)  =  0  they  are  not  ignored 

A2  =  E  upper  bound  on  truncation  error.  See  upper 
bound  Equation  (10)  Appendix  A  in  the  INT  and  INTM 
subroutine 

A3  =  M  is  value  from  which  lower  bound  E  is  calcu¬ 
lated  LBE  =  UBE/M  in  subroutine  INT 

A4  =  A  as  used  in  truncation  error  test  EO  (10)  in 
subroutine  INT 

A5  =  upper  bound  on  mesh  size  (If  =  0  no  upper  bound 
as  long  as  within  error  range) 

A6  =  lower  bound  on  mesh  size  (If  =  0  lower  bound 
=  0) 


A7  =  p,  that  is,  0  is  less  than  p  less  than  1.  It  is 
used  to  decrease  or  increase  mesh  size  by  dividing 
or  multiplying  current  integration  mesh  being  used 

smallest  attenuation  to  be  considered 

initial  refraction  index  =  W ( 1 ) 

initial  absorption  kappa  «  =  W(2) 


Table  2.  Nomenclature  Describing  the  V  and  W  Vectors. 
(Foge  4  of  6) 


51 


RM  61TMP-32 


W(51) 

W(52) 

W(53) 

W{54) 

W(55) 

W(56) 
W(57 ) 
W(58) 
W(59) 
W  (60) 
W(61) 
W(62) 
W(63) 

W  (64) 
W(65) 
W  (66) 
W(67) 

vV  (68) 


52 


initial  attenuation  =  A 

xr  (km)  (Radar  coordinate  in  geomagnetic  coordinate 
system) 

y^  (km)  (Radar  coordinate  in  geomagnetic  coordinate 
system 

(km)  (Radar  coordinate  in  geomagnetic  coordinate 
system 

R  =  J(xR  -  x)2  +  (yR  -  y)2  +  (zR  -  z)2  (km) 

Ar  =  CT  -  W (55)  =  (2.  99791  x  105)[  V(ll)  ]  -  W(55)  km 

new  elevation  angle  E  in  degrees 

Ae  =  W(57)  -  W(  17)  degrees 

2[(W(1))(W(2))]/  W(l)2  -  W(2)2 

slant  range  at  r,  9,<p 

angle  A  at  r,  9,9 

elevation  angle  E  at  r,  9,9 

assigned  value  to  k;  if  1  then  control  is  on  radius; 
if  2  then  control  is  on  range  W(28);  if  3  then  control 
is  on  slant  range  W(60) 

value  of  Z 

value  of  Y 

value  of  X 

location  of  sign  which  determines  the  calculation  for 
ordinary  or  extraordinary  ray 

value  of  V(4)  above  which  RINDEX  is  to  print  R  vector 


Table  2.  Nomenclature  Describing  the  V  and  W  Vectors. 
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W(69)  value  W(l)  below  which  R  vector  is  printed  if  W(68) 

=  0 

vV(70)  number  of  performed  integrations 

W{7  1)  a  in  COLFRZ  100-200  km 

W (72)  b  in  COLFRZ  100-200  km 

W (7 3)  a  in  COLFRZ  200-300  km 

W (74)  b  in  COLFRZ  200-300  km 

W(75)  a  in  COLFRZ  300-400  km 

W (76)  b  in  COLFRZ  300-400  km 

W(77)  to  W(250)  unassigned  in  this  program 


Table  2.  Nomenclature  Describing  the  V  and  W  Vectors. 
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B.  FUNCTION  SIANTR 


This  function  is  used  to  calculate  the  non-deviated  ray  path  between 
the  transmitter,  Rt  and  the  starting  point,  S,  (See  Figures  1  and  2) 
from  the  given  input  data  describing  the  problem  under  consideration. 
The  input  data  should  be  designed  in  a  manner  that  this  assumption 
is  true. 
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C.  FUNCTION  QATAN 


Function  QATAN  permits  the  evaluation  of  the  value  of  the  arctangent 
of  an  angle  with  proper  quadrant  allocation. 
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RM  61 TMP-32 


D.  FUNCTION  ARCOS 


Function  ARCOS  permits  the  evaluation  of  the  value  of  the  arccosine 
of  an  angle.  Presently  no  quadrant  allocation  is  made. 
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E.  SUBROUTINE  COORD 


This  subroutine  is  used  to  transform  the  problem  from  the  radar 
system  to  an  earth  centered  geomagnetic  spherical  coordinate  sys- 
tern. 
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F.  SUBROUTINE  DAUX 


Subroutine  DAUX  is  used  to  define  the  differential  equations  that 
are  to  be  numerically  integrated.  As  a  result  the  previously  des¬ 
cribed  ray  trace  equations  are  defined  in  this  subroutine. 
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G.  SUBROUTINES  INT  AND  INTM 


?; 


This  is  a  generally  available  SHARE  program  which  permits  the 
numerical  integration  of  a  chosen  set  of  first  order  non-linear  dif¬ 
ferential  equations.  It  can  be  operated  in  three  possible  numerical 
integration  modes  (a)  Runge-Kutta  with  a  fixed  integration  mesh 
size,  (b)  fourth  order  Adams -Moulton  with  a  fixed  integration  mesh 
size,  and  (c)  fourth  order  Adams-Moulton  with  a  variable  integration 
mesh  size  that  is  controlled  by  an  error  sensing  routine.  Because 
no  FORTRAN  source  program  listing  is  available,  a  SHARE  descrip¬ 
tion  of  this  routine  is  given  along  with  a  SAP  listing. 
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IDENTIFICATION 

RWINT,  Adams-Moulton,  Runge-Kutta  Integration 
704  -  FORTRAN  SAP  Language  Subroutine 
Space  Technology  Laboratories , 

Robert  Causey  and  Werner  L.  Frank, 

November  30,  1958 

ABSTRACT 

FORTRAN  version  of  RW-DE2F  which  integrates  a  system  of  N  simul¬ 
taneous,  first  order,  ordinary  differential  equations.  Option  of  using 
either  4th  order  Runge-Kutta  method  or  4th  order  predictor-corrector 
method  (Adams-Moulton)  is  provided.  Also  option  of  automatic  error 
control  with  variable  step- si  ae  is  provided.  Input  and  output  are 
single  precision  but  double  precision  is  used  internally  to  control 
round-off  errors.  Requires  12N  4  3  cells  for  data  and  693  words  for 
program. 

PURPOSE 

This  FORTRAN  subprogram  integrates  a  set  of  N  simultaneous,  first 
order  differential  equations.  It  is  the  FORTRAN  version  of  the  stand¬ 
ard  subroutine  RW-DE2F. 

RESTRICTIONS 

This  program  has  two  distinct  entries,  one  for  set  up  and  the  second 
for  performing  the  integration  steps.  The  user  must  supply  a  FOR¬ 
TRAN  subprogram  (with  the  name  DAUX)  which  evaluates  the  deriva¬ 
tives  y'  . 
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METHOD 

The  user  has  the  option  of  using  either  a  fourth  order  Runge-Kutta 
method  or  the  fourth  order  Adams-Moulton  method  with  a  fixed  step- 
size.  There  is  also  a  variable  step-size  mode. 

While  input  and  output  to  this  routine  are  single  precision,  double 
precision  is  used  internally  to  control  round-off  errors.  Truncation 
error  is  controlled  either  by  choosing  an  appropriate  step-size,  or 
by  using  the  variable  step-size  mode  of  operation. 

For  details  of  the  method  see  RW-DE2F. 

USAGE 

a.  Calling  Sequence  for  set  up  (performed  prior  to  initiating  the  in¬ 
tegration)  . 

CALL  INT  (V,  N,  Al,  A2,  A3,  A4,  A5,  A6.  A7) 


Where  V 

is  a  region  of  at  least  dimension 

N 

is  the  number  of  equations 

Al 

is  the  option  word 

A2 

is  E 

A3 

is  M 

A4 

is  A 

A5 

*8  ^max 

A6 

*8  ^min 

A7 

is  P 

For  meaning  of  Al  -  A7  see  Appendix  A  and  B  of  RW-DE2F . 

Region  V  contains  the  following  information  prior  to  Set  Up  entry. 

V  (2)  = 

V(3)  = 

V(4)  = 


V(3+N) 


x,  initial  value  of  independent  variable 
h,  value  of  step-size 

yl 


I 


=  y 


N 


values  of  dependent  variables  y 
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V(4+N)  =  y’ 
V(3+2N)  =  yN 


values  of  the  derivatives  y  to  be  supplied  by 
the  auxiliary  DAUX. 


Note:  This  region  and  the  parameter  N  should  be  placed  in 

COMMON  since  it  is  necessarily  referred  to  in  the  main 
program  and  in  the  auxiliary.  The  cell  V  (1)  is  set  up 
by  the  subprogram  RW  1NT  and  will  contain  N  scaled 
at  35. 


b.  Calling  Sequence  for  integrating  one  step. 

CALL  1NTM 

No  arguments  are  required  for  this  statement. 


SPACE  REQUIRED 

693  cells 


CHECKOUT 

This  routine  has  been  extensively  tested  on  several  check  problems. 
In  all  cases  the  errors  were  approximately  equal  to  their  expected 
values,  and  there  were  no  indications  that  round-off  errors  accumu¬ 
late  rapidly. 
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METHOD 

References: 

1.  S.  D.  Conte  and  J.  Titus,  An  interpretive  floating  point  sub¬ 
routine  for  the  solution  of  systems  of  ordinary  differential  equa¬ 
tions,  Appendix  I,  Proc.  Math.  Committee  of  Univac  Scientific 
Exchange  Meeting,  Nov.  21-22,  1957  (Obtainable  from  Reming¬ 
ton  Rand  Univac,  St.  Paul,  Minnesota). 

2.  E.  K.  Blum,  A  modification  of  the  Runge-Kutta  fourth-order 
method,  Appendix  H,  Proc.  Math.  Committee  of  Univac  Scien¬ 
tific  Exchange  Meeting,  Nov.  21-22,  1957. 


In  this  routine  the  user  is  allowed  an  option  of  using  either  the  Runge- 
Kutta  classical  fourth-order  method  as  modified  by  E.  K.  Blum 
[Ref.  (2)]  or  the  Adams -Moulton  predictor-corrector  method  using 
the  Runge-Kutta  method  for  starting  the  process.  Let  the  system  of 
equations  to  be  solved  be  given  in  the  form 


(1) 


yj  =  x-  y.> 


y{  *  =  y^ 


1  o 


io 


i  =  1,  2,  .  .  .,  N. 


Let  yin  be  the  value  of  yj  at  x  =  xn  and  fjn  the  derivation  of  y^  at 
x  =  xn  and  let  h  be  the  increment  (step-size)  of  the  independent 
variable  x.  The  classical  Runge-Kutta  fourth-order  method  uses 
the  formulas 


kil 

=  hf. 
i 

Iv 

yin| 

# 

ki2 

=  hf. 
i 

K 

t  -j  h  , 

y.  +  -rk 
in  2  ilj 

ki3 

=  hf. 

l 

1% 

+  . 

y.  +  t*.-,' 

in  2  i2 
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ki4 

=  hfi  K 

+  h  , 

yin+ki3!  ’ 

,2)  y 

. ,  +ii 

fk  ,  + 

2k  , 

+  2k,  ,  +  k  . 

yi,n+l 

yn  6  1 

1  il 

i2 

i3  14 

The  following  formulas  (we  omit  the  subscript  i  for  notational  sim¬ 
plicity)  were  derived  by  E.  K.  Blum  to  control  the  growth  of  round¬ 
off  errors. 


(3) 


(4) 


(5) 


o  n 

=  q 


4n 


P  =  h  f 


lx  *  z~ 
(  n  o 


,,  =  lO’r'1' 


Ip„  •  % 


z .  =  z  +  r 


*1 


3rl  * 


1  ’ 

—  P  -  q 
2  o  4o| 

1 


» 

> 

a 


fPl  '  h,'xnt2hl  V  ' 


r2  *  L<V2' 


z  a  z  +  r 
2  1  2 


ip.-i-.]  ■ 


t"2  *  "VT'i  *lpl 


(6) 


P2  =  2h-  *2 

N  • 


r3  •  L(3|R|3» 


*3  =  '2  3  r3 


*3  =  *  r3  +  S2 
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f  P. 


(7) 


n+ 1 


<q4,  n+1 


h  f  x  +  h,  z  I-  l  P_  , 
n  3/  Z 


l<V4> 


6P3  +  '3 


"4  "  liP3  +  q3 


,  _  (m)  (m)  ,  , 

where  R  ,  L  denote  operators  which  shift  right  m  places  or 

left  m  places  respectively  and  q^Q  is  taken  to  be  zero  to  start  the 

computation.  (See  Ref  (2)  for  a  complete  description  of  this  method.) 

Formulas  (3)  -  (7)  are  those  used  in  this  routine. 


The  Adams-Moulton  predictor-corrector  formulas  for  the  system  (1) 
are 


(8) 


y; 


(p) 


i,  n+  1 

M) 


=  y.  +  t;  55f 

in  24  1 


n'.'ntl  =  ^in  +  34 


[55  f.  -  59  f.  .  +  37  f.  ,  -  9f.  ,1  , 

|  in  i,  n-1  l,  n-2  i,  n-3J 

rz  (9f!P>  j.1  +  19  f,  -  5f.  .  +f.  , 

14  1  i,n+l  in  i,n-l  i,n-2 


The  corrector  formular  (9)  is  applied  only  once  so  that  only  two  deri¬ 
vative  evaluations  are  needed  for  each  Adams-Moulton  integration 
step.  The  starting  values  needed  in  (8)  are  obtained  using  the  Runge- 
Kutta-Blum  (RKB)  method. 


The  Adams-Moulton  method  may  be  used  either  with  a  fixed  step-size 
or  with  a  variable  step-size.  The  step-size  to  be  used  in  the  variable 
mode  is  determined  as  follows.  Let 


.<P) 


.(c) 


00) 


E  ,  =  Max 

n+1 

i 


l^i ,  n+  1  ^i ,  n+ 1 


14  D. 


>i  =  m*«  {k(:U  •  *j  • 


where  A>0.  The  user  will  specify  an  upper  bound  E  on  the  truncation 
error  estimate  E^+ ^ .  This  is  equivalent  to  specifying  the  number  of 

significant  figures  which  the  user  desires  to  preserve  locally  through¬ 
out  the  integration.  There  must  also  be  specified  a  constant  M>0 
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from  which  a  lower  bound  E  =  M"  E  is  obtained.  M  should  normally 
range  from  50  to  150.  The  interval  will  then  be  decreased,  left  as  it 
is,  or  increased  according  as  the  following  inequalities  hold: 


(Ha) 

(lib) 

(11c) 


If  E  ,  E  ,  the  interval  is  reduced  to  B  h  (0<  B  <1 
n+ 1 

If  Ex  E  ,  <  E ,  the  interval  size  is  kept  fixed. 

—  n+ 1 

If  E  <  E  for  3  successive  steps,  the  step-size  is  in- 

nr  1  —  j 

creased  to  -r  h. 

P 


Normally,  the  routine  will  take  p  =  1/2,  unless  p  is  otherwise 
specified.  The  constant  A  in  (10)  is  used  to  prevent  unnecessary 
reductions  in  |h|  whenever  Jy^  is  small.  Normally  the  routine 

will  set  A  =  1.  However,  some  other  value  for  A  may  be  specified 
by  the  user  if  he  desires  to  use  some  other  characteristic  length  for 
A.  While  the  test_based  on  (10)  will  guarantee  that  the  local  error 
does  not_exceed  E,  the  cumulative  error  will  usually  exceed  E. 
Hence,  E  should  be  chosen  small  enough  to  allow  for  an  accumula¬ 
tion  of  truncation  error.  Normally  E  should  be  in  the  range 

10'8<  E«10‘3  . 


After  an  interval  is  increased,  the  program  prevents  increasing  again 
until  6  more  points  have  been  completed.  However,  the  program  may 
decrease  the  interval  as  often  as  necessary. 


Starting  values  for  the  Adams -Moulton  formulas  are  always  obtained 
using  the  RKB  method  whenever  the  interval  size  is  changed,  just  as 
at  the  beginning  of  an  integration.  Consider  the  following  diagram  of 
the  axis  of  the  independent  variable  x 


t - ’ 


If  the  values  of  the  y.  are  computed  at  the  points  pj,  p2>  and  p3 
using  the  RKB  method  and  the  truncation  error  test  (11)  calls  for  de¬ 
creasing  |h|  at  the  point  p^,  then  the  routine  returns  to  the  point  pQ 
and  again  computes  three  new  points  with  the  RKB  method  using  the 
decreased  value  of  j li |  .  If  on  the  other  hand  (ila)  holds  at  p^  and 
the  y^  at  P3  had  been  computed  using  the  AM  method,  then  the 
routine  returns  to  the  point  P3  for  a  new  start.  If  the  inequality  in 
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(11c)  is  not  satisfied  at  pj,  but  is_  satisfied  at  p^,  p 3,  and  p^, 
then  a  new  start  is  initiated  at  with  the  increased  value  of  fh)  . 

The  user  must  provide  a  starting  value  for  h  and  he  may,  if 
desired,  specify  a  maximum  value  of  |  h|  beyond  which  the  routine 
will  not  increase  |hj  and  a  minimum  value  of  |h|  below  which  it  will 
not  decrease  |h|  .  Negative  values  of  h  may  be  supplied  for  back¬ 
ward  integration. 

Both  the  RKB  method  and  the  AM  method  incorporate  round-off  con 
trol  features.  This  is  performed  in  the  RKB  method  by  carrying 
the  q's  in  formula  (3)- (7).  Ih  the  AM  method  this  is  done  by  keep¬ 
ing  the  y^n  in  double  precision  forming  the  sums  Vin+  ^in  ix\. 
(8)  and  (9)  in  double  precision.  The  derivative  calculations  are  all 
performed  in  single  precision.  Both  procedures  have  shown  to  be 
very  effective  in  controlling  the  growth  of  round-off  errors. 
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USAGE  AND  CODING  INFORMATION  (APPENDIX  B) 


There  are  two  entries  to  this  routine.  The  first  must  be  used  once 
a t  the  beginning  to  set  up  the  routine  for  integration  of  a  given  set 
o  N  differential  equations.  The  second  entry  may  be  used  any 
number  of  times  after  the  first  to  integrate  all  y^  from  x  to  x  +  h. 
The  first  entry  has  the  following  calling  sequence. 


Loc. 

Instruction 

Al-2 

TSX  DE2F,  4 

Al-1 

PZE  T,  0,  V 

A 1 

(Binary  integer) 

A2 

(Floating-point  number) 

A3 

(Floating-point  number) 

A4 

(Floating-point  number) 

A5 

(Floating-point  number) 

A6 

(Floating-point  number) 

A7 

(Floating-point  number) 

A  7  +1 

Return 

_ Comment _ 

Setup  entry 

Parameter  word  with  addresses 
Option  word  (  =  0  or  1  or  2) 

E  I  Truncation  error  testing 
Mj  information 
A 

h  1 

maxi  jjoun<jB  on  h  ,  if  any 
"min  J 

P  -  Increase  or  decrease  h  factor 


The  eight  parameter  words  have  the  following  meaning,  (Al-1):  V 
is  the  address  of  the  first  word  of  a  block  of  12N  +  3  cells,  reserved 
by  the  user,  with  the  arrangement 


Loc. 

V 

V  +  1 

V  +  2 

V  +  3 

V  +  4 


V  +  N  +  2 


Contents 

N 

x 

h 

yf 

T2 

.  > 
yN 


Fixed  point  binary  integer,  point  at  right 
Value  of  independent  variable  in  floating  point 
Value  of  step  sise  in  floating  point 


Values  of  the  y^  in  floating  point 
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Loc.  Contents 


V  +  N  +  3 

1 

yl 

V  +  N  +  4 

1 

y2 

* 

s 

> 

• 

V  +  2N  +  2 

yN 

V  +  2N  +  3 

2 

etc. 

J 

> 

Locations  where  the  user's 
auxiliary  subroutine  must 
place  the  derivatives  yj  . 


ION  cells  of  temporary 
storage 


[  Note:  If  the  Runge-Kutta  only  option  (see  under  A1  below)  is  used, 
it  is  only  necessary  to  reserve  a  block  of  4N  +  3  cells]  . 


Before  executing  the  setup  entry,  the  user  must  have  already  placed 
the  appropriate  numbers  in  cells  V  through  V  +  N  +  2. 


The  address  V  in  the  entry  point  of  an  auxiliary  subroutine  which 
the  user  must  provide  to  evaluate  the  derivatives  y^  and  store  them 
in  cells  V  +  N  +  3  through  V  +  2N  +  2  as  shown  above.  This  auxiliary 
subroutine  is  entered  by  the  calling  sequence 


Loc.  Instruction 

A 1  -  2  TSX  V,  4 

Al-1  Return 


The  setup  entry  uses  the  auxiliary  subroutine  to  evaluate  the  deriva¬ 
tives  for  the  initial  data. 

(A1 ):  The  option  word  may  have  any  one  of  three  values  which  desig¬ 
nate  three  different  modes  of  operation  for  RWINT 

A1  =0  designates  the  predictor- corrector  variable  h  mode 
A1  =  1  designates  the  fixed  h  Runge-Kutta  only  mode 
A1  =  2  designates  the  fixed  h  predictor-corrector  mode 

For  A1  =  1  or  2,  the  contents  of  A2  through  A7  may  be  arbitrary. 
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(A2):  This  cell  contains  the  upper  hound  E>0  for  the  truncation 
error  testing  done  in  the  predictor-corrector  variable  h  mode. 
(1<T8<  E  <  10-3) 

(A3):  This  cell  contains  tha  number  M>0  from  which  the  lower 
bound  _E  is  calculated.  If  A3  =  0,  M  is  set  equal  to  100. 

(A4):  This  cell  contains  the  number  A>0  used  to  designate  a  fixed- 
point  truncation  error  test  as  described  in  Appendix  A.  If  A4  =  0, 

A  is  set  equal  to  1. 

(A5):  This  cell  may  contain  the  upper  bound  hj^^  0  for  |h|  .  If 
A5  =  0,  this  means  that  there  is  to  be  no  upper  bound  for  |h|  . 

(A 6):  This  cell  may  contain  the  lower  bound  hmjn>  0  for  jh|  .  If 
A6  =  0,  this  means  that  there  is  to  be  no  lower  bound  for  ih.|  . 

(A7):  This  cell  may  contain  the  factor  1  >  P  >  0  used  to  increase  or 
decrease  jhj  .  If  A7  =  0,  P  is  set  equal  to  1/2. 

The  integration  entry  is  quite  simple  and  has  the  calling  sequence 


Loc. 

Instruction 

Al-2 

TSX  OE2F  4  1,4  Integration  entry 

Al-1 

Return 

Ordinarily,  after  execution  of  the  integration  entry,  all  yj  assume 
new  values,  x  will  have  been  advanced  to  the  value  x  +  h  and  h 
Will  be  unchanged'  However,  in  the  variable  h  mode,  three  other 
things  can  happen.  (I)  If  the  truncation  error  test  indicates  that  jhj 

iihould  be  increased,  h  will  have  been  changed  to  P“*h  unless 
P~*hj  >hm#x.  If  the  truncation  error  test  indicates  that  |h|  should 
»e  decreased,  then  h  will  have  been  changed  to  Ph  unless 
PhJ<hm}n  and  either  (2)  y^  and  x  will  remain  as  they  were  before 
entry  or  (3)  x  will  be  changed  to  x  -  3h  and  the  corresponding  y, 
values  will  occupy  cell*  V  ♦  3  through  V  +  N  +  2.  Case  (3)  can  only 
happen  when  puecessive  decreases  in  |h|  are  called  for.  On  exit 
the  values  yj  in  V  +  N  4  3  etc.  are  always  those  which  correspond 
to  the  x  and  yj  in  V  4  1  and  V  4  3  etc. 

The  integration  entry  must  be  used  for  each  integration  step.  In 
the  variable  h  mode,  a  particular  integration  step  may  involve 
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either  AM  or  RKB  integration  but  not  both.  In  the  fixed  h  predictor- 
corrector  mode,  the  first  three  integration  entries  involve  RKB 
integration  and  all  subsequent  ones  involve  AM  integration. 

Whenever  an  integration  step  involves  AM  integration,  the  truncation 
error  estimate  En+j  is  in  the  accumulator  on  exit.  Zero  is  always 
placed  in  the  accumulator  if  the  step  involved  RKB  integration. 

The  setup  entry  may  be  used  again  at  any  time  to  set  up  another 
problem  or  to  change  the  mode  of  operation. 

In  addition  to  the  auxiliary  subroutine  for  derivative  evaluation  and 
the  12N  +  3  cells  for  data  storage,  the  storage  requirements  are 
693  words  for  RW1NT  plus  2  words  of  COMMON. 
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ORG  0 

1NT  OOOl 

DAUX 

BCD  10AUX 

INT  0002 

HTR 

|NT  0003 

HTR 

INT  0004 

HTR 

INT  0005 

1NT 

SXD  INT-3.4 

INT  0006 

SXD  I  NT-2 *2 

INT  0007 

SXD  INT-1.1 

INT  0008 

CLA  1 *4 

INT  0009 

STA  REV1+1 

INT  0010 

STA  REVl+2 

INT  0011 

STA  A2 

INT  0012 

CLA  2*4 

INT  0013 

STA  A1 

INT  0014 

A1 

CLA  N 

INT  0015 

ARS  18 

INT  0016 

A2 

STO 

INT  0017 

ALS  2 

INT  0018 

STO  C 

INT  0019 

ALS  1 

INT  0020 

ADD  C  -12M 

INT  0021 

ADD  Cl  12N+2 

INT  0022 

STO  C 

INT  0023 

CLA  1.4 

INT  0024 

SUB  C 

INT  0023 

STA  PARI 

INT  0026 

CLA  C 

INT  0027 

ARS  1 

INT  0028 

STO  C 

INT  0029 

ADD  PARI 

INT  0030 

STA  REV1 

INT  0031 

STA  REVl+3 

INT  0032 

LXA  C2  *1 

INT  0033 

A4 

CLA  3.4 

INT  0034 

STA  AS 

INT  0035 

AS 

CLA 

INT  0036 

STO  PAR8+1 *1 

INT  0037 

TXI  A6.4.-1 

INT  0038 

A6 

TIX  A4.1.1 

INT  0039 

CLA  PAR2 

INT  0040 

LRS  18 

INT  0041 

STO  PAR2 

INT  0042 

CLA  OAUX 

INT  0043 

STA  AUX+2 

INT  0044 

TSX  REV. 4 

INT  0045 

TSX  DE2P  *4 

INT  0046 

PARI 

PZE  O.O.AUX 

INT  0047 

PAR  2 

P2E 

INT  0048 

PARS 

PZE 

INT  0049 

PARA 

PZE 

INT  0050 

PARS 

PZE 

INT  0051 

PARS 

PZE 

INT  0052 

PART 

PZE 

INT  0053 
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PARS 

P2E 

I  NT 

0054 

TSX 

REV.4 

INT 

0055 

LXD 

I  NT-2  »2 

INT 

0056 

LXD 

INT-1  »  1. 

INT 

0057 

LXD 

I NT-3  *4 

INT 

0058 

TRA 

10«4 

INT 

0059 

1NTM 

SXD 

INT-3.4 

INT 

0060 

SXD 

I  NT-2  *2 

INT 

0061 

SXD 

1  NT-1  *  1 

INT 

0062 

TSX 

REV. 4 

INT 

0063 

TSX 

DE2F+1 «4 

INT 

0064 

TSX 

REV.4 

INT 

0065 

LXD 

INT-3.4. 

INT 

0066 

LXD 

INT-2.2 

INT 

0067 

LXD 

INT-1. 1 

INT 

0068 

TRA 

1.4 

INT 

0069 

REV 

LXA 

C.l 

INT 

0070 

LXD 

Cl. 2 

INT 

0071 

REV1 

CLA 

0.1 

INT 

0072 

LDO 

0.2 

INT 

0073 

STO 

0.2 

INT 

0074 

STO 

0.1 

INT 

0075 

TXI 

REV3.2.1 

INT 

0076 

REV3 

TIX 

REV1.1.1 

INT 

0077 

TRA 

1  .4 

INT 

0078 

AUX 

SXO 

C3.4 

INT 

0079 

TSX 

REV. 4 

INT 

0080 

TSX 

0.4 

INT 

0081 

TSX 

REV. 4 

INT 

0082 

LXD 

C3.4 

INT 

0083 

TRA 

1.4 

INT 

0084 

C 

PZE 

INT 

0085 

Cl 

DEC 

2 

INT 

0086 

C2 

OEC 

7 

INT 

0087 

C3 

INT 

0088 

REM 

FLOATING  POINT 

ADAMSr-HOULTON .  RUNGE-KUTTA  INTEGRATION 

INT 

0089 

DE2F 

TRA 

DE2F+0293 

SETUP  ENTRY 

INT 

0090 

SXD 

DE2F+0240»1 

INTEGRATION  ENTRY 

INT 

0091 

SXD 

DE2F+0241.2 

INT 

0092 

SXD 

DE2F+0242 .4 

INT 

0093 

TRA 

DE2F+0005 

SWITCH  1 

INT 

0094 

CLA 

DE2F+0263 

3  TO  ACC 

INT 

0095 

CAS 

0E2F+022S 

TEST  ALPHA 

INT 

0096 

TRA 

DE2F+0186 

INT 

0097 

TRA 

DE2F+0175 

INT 

0098 

LXA 

DE2F+0229. 1 

Y  PRIMES  TO  D 

INT 

0099 

CLA 

*.l 

INT 

0100 

STO 

•  .1 

AND 

INT 

0101 

CLA 

•  .1 

INT 

0102 

STO 

».l 

DOUBLE  PRECISION 

INT 

0103 

CLA 

*.l 

INT 

0104 

STO 

*.l 

YS  TO  TS2 

INT 

0105 

TIX 

DE2F+0010.1.1 

END  3F  LOOP 

INT 

0106 
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LXA  DE2F+0229 tl 
LOO  DE2F+0223 
FMP  +  *1 

STO  COMMON+OOO 

CLA  +  .1 

STO  DE2F+0267 

CLA  *tl 

STO  DE2F+0268 

LOO  DE2F+0222 

FMP  +  «1 

FAD  COMMON+OOO 
STO  COMMON+OOO 
LOO  OE2F+0221 
FMP  +  .1 

FAD  COMMON+OOO 
STO  COMMON+OOO 
LDO  DE2F+0220 
FMP  *.l 

FAD  COMMON+OOO 
STO  COMMON+OOO 
LDO  * 

FMP  COMMON+OOO 
TSX  DE2F+0269.2 
STO  *.l 

T1X  DE  2F+0018  *  1*1 
CLA  * 

FAD  * 

STO  • 

LXA  0E2F+0229.1 
CLA  *tl 
STO  ».l 

T1X  DE2F+0046tl»l 
TSX  0.4 

LXA  DE2F+0229#! 
LOO  DE2F+0227 
FMP  *.l 

STO  COMMON+OOO 

CLA  *.l 

STO  DE2F+0267 

CLA  +  .1 

STO  OE2F+0268 

LDO  DE2F+0226 

FMP  •  »! 

FAD  COMMON+OOO 
STO  COMMON+OOO 
LOO  DE2F+0225 
FMP  **1 

FAD  COMMON+OOO 
STO  COMMON+OOO 
LOO  DE2F+0224 
FMP  »*1 

FAD  COMMON+OOO 
STO  COMMON+OOO 


D3  SUB  I 

PLANT  Y  SUB  1 

IN  SFA 

SUBROUTINE' 

02  SUB  t 

01  SUB  I 

D  SUB  ! 

LOAD-  H 

DELTA  YI  UPPER  ► 

AOD  TO  YI 
STORE  IN  TS1 
END  OF  LOOP 

X+H 

ENO  OF  LOOP 
EVALUATE  DERIVATIVES 

D2  SUB  I 
YI  FROM  TS2 
DP  EXT* 

D1  SUB I 

d  suet 

Y  PRIME  SUB! 


INT  0107 
INT  0108 
INT  0109 
INT  0110 
INT  0111 
INT  0112 
INT  0113 
INT  0114 
INT  0115 
INT  0116 
INT  0117 
INT  0118 
INT  0119 
INT  0120 
INT  0121 
INT  0122 
INT  0123 
INT  0124 
INT  0125 
INT  0126 
INT  0127 
INT  0128 
INT  0129 
INT  0130 
INT  0131 
INT  0132 
INT  0133 
INT  0134 
INT  0135 
INT  0136 
INT  0137 
INT  0138 
INT  0139 
INT  0140 
INT  0141 
INT  0142 
INT  0143 
INT  0144 
INT  0145 
INT  0146 
INT  0147 
INT  0148 
INT  0149* 
INT  0150 
INT  0151 
INT  0152 
INT  0153 
INT  0154 
INT  0155 
INT  0156 
INT  0157 
INT  0158 
INT  0159 
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LOO 

* 

LOAD  H 

INT 

0160 

FMP 

COMMON+OOO 

DELTA  YI  UPPER  C 

INT 

C 1 6 1 

TSX 

DE2F+Q269  *  2 

ADD  TO  Y  I 

INT 

0162 

STO 

*.  1 

INT 

0163 

STO 

*.l 

INT 

0164 

T  I X 

DE2F+0051 *1*1 

END  OF  LOOP 

INT 

0165 

STZ 

DE2F+0230 

INT 

0166 

LXA 

DE2F+0229 • 1 

INT 

0167 

CLA 

*,1 

INT 

0168 

STO 

DE2F+0243 

INT 

0169 

CLA 

**  1 

INT 

0170 

STO 

DE2F+0244 

INT 

0171 

TSX 

DE2F+0246.2 

INT 

0172 

T I X 

DE2F  +  0078  *1*1 

END  OF  LOOP 

INT 

0173 

TRA 

DE2F+0085 

SWITCH  2 

INT 

0174 

CLA 

DE2F+0234 

INT 

0175 

CAS 

DE2F+0230 

INT 

0176 

TRA 

DE2F+0090 

INT 

0177 

TRA 

DE2F+0089 

INT 

0178 

TRA 

DE2F+0115 

DECREASE  H  SWITCH 

INT 

0179 

CLA 

DE2F+0235 

INT 

0180 

CAS 

DE2F+0230 

INT 

0181 

TRA 

OE2F+0164 

INCREASE  H  SWITCH 

INT 

0182 

TRA 

DE2F+0211 

INT 

0183 

TRA 

DE2F+02 1 1 

INT 

0184 

LXA 

DE2F+0229.1 

INT 

0185 

CLA 

*.l 

INT 

0186 

STO 

*.l 

D2  TO  03 

INT 

0187 

CLA 

*.l 

INT 

0188 

STO 

*.l 

D1  TO  D2 

INT 

0189 

CLA 

*.l 

INT 

0190 

STO 

•*1 

D  TO  D1 

INT 

0191 

T  I X 

OE  2F+0096  *1*1 

END  OF  LOOP 

INT 

0192 

CLA 

DE2F+0228 

INT 

0193 

ADD 

0E2F^0284 

ALPHA  PLUS  ONE 

INT 

0194 

STO 

0E2F+0228 

INT 

0195 

TSX 

0*4 

EVALUATE  DERIVATIVES 

INT 

0196 

CLA 

0E2F+0230 

INT 

0197 

FOP 

0E2F+0290 

E 

INT 

0198 

STO 

COMMON+OOO 

INT 

0199 

CLA 

COMMON+OOO 

GET  E  IN  ACC 

INT 

0200 

LXO 

0E2F+0240* 1 

INT 

0201 

LXO 

OE2F+0241.2 

INT 

0202 

LXO 

0E2F+0242  *4 

INT 

0203 

TRA 

1*4 

EXIT 

INT 

0204 

CLA 

• 

INT 

0203 

SSP 

INT 

0206 

CAS 

OE2F+0238 

TEST  H  WITH  HMIN 

INT 

0207 

TRA 

DE2F+0121 

INT 

0208 

TRA 

OE2F+0121 

INT 

0209 

TRA 

0E2F+0211 

INT 

0210 

LOO 

• 

INT 

0211 

STO 

DE2F+0233 

STORE  OLD  H 

INT 

0212 
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FMP  DE2F+0239 

BETA  TIMES  H 

1NT  0213 

STO  * 

INT  0214 

CLA  OE2F+0285 

INT  0215 

CAS  OE 2F+0228 

TEST  ALPHA 

INT  0216 

HTR  OE 2F+0127 

INT  0217 

TRA  DE2F+0130 

INT  0218 

TRA  DE2F+0144 

INT  0219 

LOO  DE2F+0289 

INT  0220 

FMP  DE2F+0233 

4H 

INT  0221 

CHS 

INT  0222 

FAO  * 

X-4H 

INT  0223 

STO  * 

INT  0224 

LXA  DE2F+0229.1 

INT  0225 

CLA  *.l 

INT  0226 

STO  4.1 

INT  0227 

CLA  4.1 

INT  0228 

STO  *.l 

INT  0229 

CLA  *.l 

INT  0230 

STO  *.l 

INT  0231 

TlX  DE2F+0136.1.1 

END  OF  LOOP 

INT  0232 

TRA  DE2F+0161 

JUMP  TO  SET  ALPHA 

INT  0233 

CLA  DE2F+0233 

INT  0234 

CHS 

INT  0235 

FAD  * 

X-H 

INT  0236 

STO  * 

INT  0237 

LXA  OE2F+0229.1 

INT  0238 

CLA  *. 1 

INT  0239 

STO  *.l 

0  TO  YI  PRIME 

INT  0240 

CLA  *.l 

T52  TO  YI 

INT  0241 

STO  *.l 

INT  0242 

CLA  *.l 

INT  0243 

STO  *.l 

INT  0244 

TlX  DE2F+0149.1.1 

INT  0245 

LXA  DE2F+0229.1 

INT  0246 

LOO  OE2F+0291 

CONVERT  OP  YI 

INT  0247 

FMP  *.l 

TO  YI.  01 

INT  0248 

STO  *.l 

INT  0249 

T  t X  OE2F+0157.1.1 

END  OF  LOOP 

INT  0250 

ST2  OE2F40228 

INT  0251 

STZ  0E2F40231 

INT  0252 

TRA  DE2F+0107 

INT  0253 

CLA  * 

INT  0254 

SSP 

INT  0255 

CAS  DE2F+0237 

TEST  H  WITH  HMAX 

INT  0256 

TRA  DE2F+0095 

GO  TO  SHIFT  DS 

INT  0257 

TRA  0E2F40213 

INT  0250 

TRA  OE2F40213 

INT  0259 

CLA  * 

INT  0260 

FOP  DE2F+0239 

H  DIVIDED  BY  BETA 

INT  0261 

STO  * 

INT  0262 

tsx  O.A 

evaluate  derivatives 

INT  0263 

TRA  DE2F+6156 

M  °?6t 

LxA  DE2F*0229il 

INT  0265 
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CUA  *,1 

INT  0266. 

FOP  DE2F+0291 

0!  DIVIDED  BY  -3 

INT  0267 

CIA  *,1 

INT  0268 

STO  COMMON+OOO 

INT  0269 

FAO  COMMON+OOO 

INT  0270 

STO  *.l 

INT  0271 

STO  *«1 

INT  0272 

TIX  DE2F+0176.1.1 

INT  0273 

TRA  DE2F+0009 

INT  0274 

STO  *.l 

DUMMY  ORDER 

INT  0275 

LOQ  DE2F+0229 

INT  0276 

MPY  DE2F+0228 

INT  0277 

STQ  COMMON+OOO 

INT  0278 

CLA  COMMON+OOO 

INT  0279 

ADD  OE2F+0185 

INT  0280 

STO  DE2F+0194 

INT  0281 

LXA  0E2F+0229 *  1 

INT  0282 

CLA  *.l 

INT  0283 

STO  0 

STORE  DERIVATIVES 

INT  0284 

TIX  DE2F+0193.1.1 

INT  0285 

CLA  OE2F+0228 

INT  0286 

TNZ  DE2F+0204 

INT  0287 

LXA  DE2F+0229  tl 

INT  0288 

CLA  *.l 

INT  0289 

STO  *.l 

INT  0290 

CLA  *.l 

INT  0291 

STO  *.l 

INT  0292 

TIX  DE2F+0199.1.1 

INT  0293 

CLA  OE2F+0228 

INT  0294 

ADO  DE2F+0284 

ALPHA  PLUS  ONE 

INT  0295 

STO  DE2F+0228 

INT  0296 

TSX  DE2F+0454  *4 

RUNGA-ICUTTA  ENTRY 

INT  0297 

STZ  COMMON+OOO 

INT  0298 

CLA  COMMON+OOO 

ZERO  TO  ACC 

INT  0299 

TRA  DE2F+0111 

INT  0300 

STZ  DE2F+0231 

INT  0301 

TRA  DE2F+0095 

INT  0302 

CLA  DE2F+0231 

INT  0303 

ADD  0E2F+0284 

R+l 

INT  0304 

STO  0E2F+0231 

INT  0305 

CAS  0E2F+0292 

INT  0306 

HTR  DE2F+0217 

INT  0307 

TRA  DE2F+0170 

INT  0308 

TRA  DE2F+0095 

INT  0309 

DEC  2.291666667 

INT  0310 

DEC  -2.458333333 

INT  0311 

OEC  1.541666667 

INT  0312 

DEC  -3.75E-1 

INT  0313 

DEC  3.75E-1 

INT  0314 

DEC  7.916666667E-1 

INT  0315 

OEC  -2.083333333E-1 

INT  0316 

DEC  4. 166666667E-2 

INT  0317 

BSS  4 

INT  0318 
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BSS 

8 

I  NT 

0319 

ess 

3 

INT 

0320 

HTR 

HI  UPPER  C 

INT 

0321 

HTR 

YI  UPPER  P 

INT 

0322 

HTR 

DIVISOR 

INT 

0323 

CLA 

DE2F+0243 

ENTRY 

INT 

0324 

SSP 

INT 

0323 

CAS 

DE2F+0236 

TEST  FOR  DIVISOR 

INT 

0326 

TRA 

DE2F+0254 

INT 

0327 

TRA 

DE2F+0254 

INT 

0328 

CLA 

DE2F+0236 

INT 

0329 

STO 

DE2F+0245 

INT 

0330 

TRA 

DE2F+0255 

INT 

0331 

STO 

0E2F+0245 

INT 

0332 

CLA 

DE2F+0244 

INT 

0333 

FSB 

DE2F+0243 

INT 

0334 

SSP 

INT 

0335 

FDP 

DE2F+0245 

INT 

0336 

STO 

COMMON+OOO 

INT 

0337 

CLA 

DE2F+0230 

INT 

0338 

CAS 

COMMON+OOO 

INT 

0339 

TRA 

1.2 

INT 

0340 

TRA 

1.2 

INT 

0341 

CLA 

COMMON+OOO 

INT 

0342 

STO 

DE2F+0230 

INT 

0343 

TRA 

1.2 

INT 

0344 

HTR 

A1 

INT 

0345 

HTR 

A2 

INT 

0346 

UFA 

DE2F+0267 

ENTRY 

INT 

0347 

STO 

0E2F+0267 

INT 

0348 

STO 

COMMON+OOO 

SPECIAL 

INT 

0349 

CLA 

COMMON+OOO 

FLOATING 

INT 

0350 

UFA 

DE2F+0268 

ADDITION 

INT 

0351 

FAD 

0E2F+0267 

SUBROUTINE 

INT 

0352 

TRA 

1.2 

INT 

0353 

TRA 

DE2F+0005 

SWITCH  1*  A  LEG 

INT 

0354 

TRA 

DE2F+0207 

SWITCH  1.  B  LEG 

INT 

0355 

TRA 

OE2F+0O83 

SWITCH  2.  A  LEG 

INT 

0356 

TRA 

DE2F+0095 

SWITCH  2.  B  LEG 

INT 

0357 

TRA 

0E2F+0115 

DECREASE  H  SWITCH* 

TEST 

LEG 

INT 

0358 

TRA 

OE2F+0121 

DECREASE  H  SWITCH* 

NO  TEST  LEG 

INT 

0359 

TRA 

0E2F+0164 

INCREASE  H  SWITCH* 

TEST 

LEG 

INT 

0360 

TRA 

0E2F+0213 

INCREASE  H  SWITCH* 

NO  TEST  LEO 

INT 

0361 

DEC 

1 

INT 

0362 

OEC 

3 

INT 

0363 

OEC 

1. 

INT 

0364 

OEC 

5E-1 

ONE  HALF 

INT 

0365 

OEC 

100. 

INT 

0366 

OEC 

4. 

INT 

0367 

OEC 

14. 

INT 

0368 

OEC 

-3. 

INT 

0369 

DEC 

3 

SPECIAL  TEST  NO.  FOR  R 

INT 

0370 

SXD 

DE2F+0241.2 

INT 

0371 
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SXD  DE2F+0242  *4 
CLA  DE2F+0276 
STO  DE2F+0004 
CLA  DE2F+0278 
STO  DE2F+0084 
CLA  0E2F+0280 
STO  DE2F+0089 
CLA  0E2F+0282 
STO  DE2F+0092 
CLA  2.4 
PAX  0.2 

TRA  DE2F+0308 . 2 
TRA  DE2F+0346 
TRA  DE2F+0343 
LOO  3.4 
FMP  DE2F+0290 
STO  0E2F+0234 
CLA  8.4 
TNZ  DE2F+0314 
CLA  DE2F+0287 
STO  DE2F+0239 
CLA  4,4 
TNZ  DE2F+03 18 
CLA  DE2F+0288 
STO  DE2F+0235 
CLA  DE2F+0234 
FOP  DE2F+0235 
STO  DE2F+0235 
CLA  5.4 
TNZ  0E2F+0325 
CLA  DE2F+0286 
STO  OE2F+0236 
CLA  6.4 
TNZ  0E2F+0331 
CLA  DE2F+0283 
STO  OE2F+0092 
TRA  DE2F+0335 
STO  DE2F+0237 
LOO  OE2F+0239 
FMP  DE2F+0237 
STO  DE2F+0237 
CLA  7,4 
TNZ  0E2F-.0340 
CLA  DE2F+0281 
STO  DE2F+0089 
TRA  DE2F+0348 
FOP  DE2F*0239 
STO  DE2F+0238 
TRA  OE  2F+0348 
CLA  OE2F+0277 
STO  DE2F+0004 
TRA  OE 2F*0348 
CLA  DE2F+0279 


INT 

0372 

INT 

0373 

SET 

INT 

0374 

SWITCHES 

INT 

0375 

TO 

INT 

0376 

NORMAL 

INT 

0377 

POSITIONS 

INT 

0378 

INT 

0379 

INT 

0380 

INT 

0381 

INT 

0382 

SELECT  OPTION 

INT 

0383 

OPTION  2 

INT 

0384 

OPTION  1 

INT 

0385 

OPTION  0 

INT 

0386 

14E  (UPPER) 

INT 

0387 

INT 

0388 

INT 

0389 

TEST  BETA 

INT 

0390 

INT 

0391 

STORE  BETA 

INT 

0392 

INT 

0393 

INT 

0394 

INT 

0395 

STORE  H 

INT 

0396 

INT 

0397 

INT 

0398 

STORE  14E  (LOWER) 

INT 

0399 

INT 

0400 

TEST  A 

INT 

0401 

INT 

0402 

STORE  A 

INT 

f  fclT 

0403 

A/.n/. 

TEST « H  MAX 


STORE  BETA  (H  MAXI 
TEST  H  MIN 


INT  0405 
I  NT  0406 
INT  0407 
INT  0408 
INT  0409 
INT  0410 
INT  0411 
INT  0412 
INT  0413 
INT  0414 
INT  0415 
INT  0416 
INT  0417 
INT  0418 
INT  0419 
INT  0420 
INT  0421 
INT  0422 
INT  0423 
INT  0424 
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STO  DE2F+0084 
CLA  1 »4 
STO  DE2F+0232 
STO  DE2F+0352 
TSX  DE  2F  +  0453  »4 
PZE 

STZ  DE2F+0228 
STZ  DE2F+0231 
CLA  DE2F+0232 
STA  DE2F+0357 
CLA  0 

STO  DE 2F+0229 
CLA  DE2F+0232 
ARS  18 

STA  DE2F+0049 
STA  DE2F+0106 
STA  0E2F+0173 
CLA  DE2F+0232 
ADO  0E2F+02B4 
STA  DE2F+0042 
STA  DE2F+0044 
STA  DE2F+0133 
STA  DE2F+0134 
STA  DE2F+0146 
STA  DE2F+0147 
ADD  OE 2F+0284 
STA  DE2F+0037 
STA  DE2F+0043 
STA  DE  2F+0070 
STA  DE2F+0115 
STA  DE2F+0121 
STA  DE2F+0124 
STA  DE2F+0164 
STA  DE2F+0170 
STA  DE2F+0172 
ADO  DE2F+02B4 
ADO  OE2F *0229 
STA  DE2F+0010 
STA  DE2F+0021 
STA  DE2F+0047 
STA  OE2F+0073 
STA  DE2F+0078 
STA  DE2F+0139 
STA  0E2F«-0152 
STA  DE2F+0178 
STA  DE2F+0181 
STA  DE2F+0199 
ADD  DE2F+0229 
STA  DE2F40012 
STA  0E2F+0067 
STA  DE2F*0137 
STA  DE2F+0150 
STA  DE2F*0193 


SETUP  RK  SUBROUTINE 
PAR 'METER  WORD 
SET  ALPHA  TO  ZERO 
R  0 


STORE  N 


SETUP 

DERIVATIVE 

EVALUATIONS 

T-l 


T-2 


T-3  EOUALS  0 

D*N 


D+2N 


INT  0425 
1NT  0426 
INT  0427 
INT  0428 
INT  0429 
INT  0430 
INT  0431 
INT  0432 
INT  0433 
INT  0434 
INT  0435 
INT  0436 
INT  0437 
INT  0438 
INT  0439 
INT  0440 
INT  0441 
INT  0442 
INT  0443 
INT  0444 
INT  0445 
INT  0446 
INT  0447 
INT  0448 
INT  0449 
INT  0450 
INT  0451 
INT  0452 
INT  0453 
INT  0454 
INT  0455 
INT  0456 
INT  0457 
INT  0458 
INT  0459 
INT  0460 
INT  0461 
INT  0462 
INT  0463 
INT  0464 
INT  0465 
INT  0466 
INT  0467 
INT  0468 
INT  0469 
INT  0470 
INT  0471 
INT  0472 
INT  0473 
INT  0474 
INT  0475 
INT  0476 
INT  0477 
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ADO  DE2F+0229 

D+3H 

INT  0478 

STA  DE2F+0014 

INT  0479 

STA  DE2F+0023 

INT  0480 

STA  DE2F+0056 

INT  0481 

STA  DE2F+0074 

INT  0482 

STA  DE2F+0141 

INT  0483 

STA  0E2F  +  0 1 54 

INT  0484 

STA  DE2F+0158 

INT  0485 

STA  DE2F+0 1 59 

INT  0486 

STA  0E2F+0176 

INT  0487 

STA  DE2F+0182 

INT  0488 

STA  DE2F+0201 

INT  0489 

AOD  DE2F+0229 

D+4N 

INT  0490 

STA  DE2F+0040 

INT  0491 

STA  DE2F+0046 

INT  0492 

STA  DE2F+0080 

INT  0493 

ADD  DE2F+0229 

D+5N 

INT  0494 

STA  DE2F+0011 

INT  0495 

STA  DE2F+0054 

INT  0496 

STA  DE2F+01 51 

INT  0497 

ADD  DE2F+0229 

D+6N 

INT  0498 

STA  DE2F+O015 

INT  0499 

STA  DE2F+0153 

INT  0500 

AOD  DE2F+0229 

0+7N 

INT  0501 

STA  0E2F+0138 

INT  0502 

STA  DE2F+0200 

INT  0503 

ADO  DE2F+02 29 

D+8N 

INT  0504 

STA  DE2F+0140 

INT  0505 

STA  0E2F+0202 

INT  0506 

ADO  DE2F+0229 

D+9N 

INT  0507 

STA  DE2F+0019 

INT  0508 

STA  DE2F+0097 

INT  0509 

STA  0E2F+Q1 36 

INT  0510 

STA  DE2F+0185 

INT  0511 

ADD  DE2F+0229 

D+10N 

INT  0512 

STA  DE2F+0026 

INT  0513 

STA  DE2F+O052 

INT  0514 

STA  DE2F+O096 

INT  0515 

STA  DE2F+0099 

INT  0516 

AOD  DE2F+0229 

D*11N 

INT  0517 

STA  OE2F+O030 

INT  0518 

STA  DE2F+O059 

INT  0519 

STA  DE2F+0098 

INT  0520 

STA  DE2F+0101 

INT  0521 

ADD  DE2F*0229 

0+12N 

INT  0522 

STA  DE2F+0013 

INT  0523 

STA  DE2F+0034 

INT  0524 

STA  DE2F+0063 

INT  0525 

STA  0E2F+0100 

INT  0526 

STA  DE2F+0149 

INT  0527 

IXD  DE2Ff024l.2 

INT  0528 

LXD  0E2F+0242.4 

INT  0529 

TRA  9.4 

EXIT 

INT  0530 
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TRA  DE2F+0562 
SXO  DE2F+0553.1 
SXD  DE  2F  +  0554  •  2 
SXD  DE2F+0555  *4 
CLA 

FDP  DE2F+0520 
STQ  DE2F+0525 
STQ  DE2F+0541 
LXA  DE2F+0551*2 
LXA  DE2F+0556* 1 
IDO  *1 

FMP  DE2F+0551 »2 

STO  .1 

LOO 

FMP  *1 
FAO  .1 
STO  .1 

FDP  DE2F+0552»2 
STO  COMMON+OOO 
LDO  .1 

FMP  DE2F+0553*  2 
FAD  COMMON+OOO 
STO  COMMON+OOO 
TZE  DE2F+0492 
ARS  27 

STO  COMMON+OOl 
CLA  .1 

TZE  DE2F+0492 
SSP 

ARS  27 

SBM  COMMON+OOl 
TM!  DE2F+0492 
TZE  DE2F+0492 
STA  DE2F+0489 
STA  0E2F+0490 
CLA  COMMON+OOO 
ARS 
ALS 

STO  COMMON+OOO 
CLA  .1 

FAO  COMMON+OOO 
STO  »1 

TRA  DE2F+0558.2 
LOO  «1 

FMP  DE2F+0556*2 
STO  tl 
LOO  .1 

FMP  DE2F+0555.2 
FAO  »1 
STO  »1 

LOO  COMMON+OOO 
FMP  DE2F+0554.2 
FAO  *1 


TO  SETUP  REGION 
SAVE  INDEX 
FROM 

MAIN  PROGRAM 


REGISTERS 


CALCULATE 
H  DIVIDED  BY  2 

SET  PARAMETER  INDEX 
SET  N  INDEX. 


CALCULATE  NEW 
VALUE  OF  P 


CALCULATE  NEW 
VALUE  OF  R 


TEST  VALUES 
TO 

DETERMINE 
IF  SHIFTING 
IS  NECESSARY 


CALCULATE  NEW 
VALUE  OF  Z 

CALCULATE 
NEW  VALUE 

OF  0 


INT  OS31 
INT  0532 
INT  0533 
INT  0534 
INT  0535 
INT  0536 
INT  0537 
INT  0538 
INT  0539 
INT  0540 
INT  0541 
INT  0542 
INT  0543 
INT  0544 
INT  0545 
INT  0546 
INT  0547 
INT  0540 
INT  0549 
INT  0550 
INT  0551 
INT  0552 
INT  0553 
INT  0554 
INT  0555 
INT  0556 
INT  0557 
INT  0550 
INT  0559 
INT  0560 
INT  0561 
INT  0562 
INT  0563 
INT  0564 
INT  0565 
INT  0566 
INT  0567 
INT  0560 
INT  0569 
INT  0570 
INT  0571 
INT  0572 
INT  0573 
INT  0574 
INT  0575 
INT  0576 
INT  0577 
INT  0570 
INT  0579 
INT  0360 
INT  0581 
INT  0582 
INT  0503 
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STO  *1 

T 1 X  DE2F+0463.1»1  TEST  N 

CLA 

FAD  DE2F+0557.2  INCREASE  X 

STO 

SXD  DE2F+0557.2 

TSX  .4  FIND  DERt VlTIVES 

LXD  DE2F*0537»2 

T IX  DE2P+0462  *2*8  TEST  PASS  NO. 

LXD  DE2P+0553 » 1  RESTORE 

LXD  DE2F+0554.2  INDEX 

LXD  DE2F+0553  *4  REGISTERS 

TRA  1,4 
DEC  0 
DEC  2. 

DEC  -1. 

DEC  3. 

DEC  -.5 
DEC  1. 

DEC 

TRA  DE2F+0496 
DEC  0 
DEC  2. 

DEC  -.5 
DEC  -1. 

DEC  .5 
DEC  -3. 

DEC  0 

TRA  DE2F+0358 
DEC  -.5 
DEC  1. 

DEC  0 
DEC  -1. 

DEC  0 
DEC  1. 

DEC 

TRA  DE2F40496 
DEC  2. 

OEC  6. 

DEC  1. 

OEC  3. 

OEC  -.3 
DEC  -3. 

OEC  0 

TRA  DE2F+0496 
P2E  32,, 

PZE  1,, 

BSS  5 
CLA  ,1 

FDP  DE2F+05S6.2 
STO  ,1 

TRA  0E2F40499 

SXO  0E2F+0553, 1  SAVE  INDEX  REGISTERS 


INT  0584 
INT  0585 
INT  0586 
INT  0587 
INT  0588 
INT  0589 
INT  0590 
INT  0591 
INT  0592 
INT  0593 
INT  0594 
INT  0595 
INT  0596 
INT  0597 
INT  0598 
INT  0599 
INT  0600 
INT  0601 
INT  0602 
INT  0603 
INT  0604 
INT  0605 
INT  0606 
INT  0607 
INT  0608 
INT  0609 
INT  0610 
INT  0611 
INT  0612 
INT  0613 
INT  0614 
INT  0615 
INT  0616 
INT  0617 
INT  0618 
INT  0619 
INT  0620 
INT  0621 
INT  0622 
INT  0623 
INT  0624 
INT  0625 
INT  0626 
INT  0627 
INT  0628 
INT  0629 
INT  0630 
INT  0631 
INT  0632 
INT  0633 
INT  0634 
INT  0635 
INT  0636 


T 


110 


RM  61 TM  P-32 


SXD  DE2F+0555.4 
CLA  1.4 
STA  DE2F+0553 
STA  DE2F+0570 
ARS  16 

STA  OE2F+0512 
STA  DE2F+0603 
CLA 

STO  DE2F+0556 
CLA  DE2F+0553 
ADD  DE2F+0552 
STA  DE2F+0508 
STA  OE2F+0510 
ADO  DE2F+0552 
STA  DE2F+0466 
STA  DE2F+0457 
ADO  DE2F+0552 
ADO  DE2F+0556 
STA  0E2F+0479 
STA  DE2F+0492 
STA  DE2F+0494 
ADD  DE2F+0556 
STA  DE2F+0467 
ADD  0E2F+0536 
STA  DE2F+0603 
STA  DE2F+0560 
STA  DE2F+0558 
STA  DE2F+047Z 
STA  DE2F+0496 
STA  DE2F+0498 
STA  0E2F+0301 
STA  DE2F+0502 
STA  DE2F-*0303 
STA  DE2F+0506 
ADD  DE2F+0536 
STA  0E2F+0463 
STA  DE2F+0463 
STA  OE2F+0468 
STA  DE2F+0469 
STA  DE2P+0499 
TSX  .4 

LXA  DE2F+0556»1 
ST2  <1 

TIX  DE2F+0605*1«1 
LXD  -E2F+05S3.1 
LXD  DE2F+0553»* 
TRA  2*4 
COMMON  BSS  2 
R  END 


FROM  MAIN  PROGRAM 


SET  ADDRESS  OF 
DERIVITlVE  ROUTINE 
STORE  VALUE 
OF  N 


STORE  ADDRESS 
OF  X 

STORE  ADDRESS 
OF  H 


STORE  ADORE SS 
OF  Y 

STORE  ADDRESS 
OF  DERIVITlVE 


STORE 

ADDRESS 

OF 

0 


STORE  ADORESS 
OF  P 

FIND. INITIAL  DCRIVITIVCS 

SET  OROINAL  0 
TO  ZERO 
RESTORE  INDEX 
REGISTERS 


I N  T  0637 
INT  0638 
INT  0639 
INT  0640 
INT  0641 
INT  0642 
INT  0643 
INT  0644 
INT  0643 
INT  0646 
INT  0647 
INT  0648 
INT  0649 
INT  0650 
INT  0651 
INT  0652 
INT  0653 
INT  0654 
INT  0655 
INT  0656 
INT  0637 
INT  0658 
INT  0659 
INT  0660 
INT  0661 
INT  0662 
INT  0663 
INT  0664 
INT  0665 
INT  0666 
INT  0667 
INT  0668 
INT  0669 
INT  0670 
INT  0671 
INT  0672 
INT  0673 
INT  0674 
INT  0673 
INT  0676 
INT  0677 
INT  0678 
INT  0679 
INT  0680 
INT  0681 
INT  0682 
INT  0683 
INT  0684 
INT  0685 
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H.  SUBROUTINE  RINDEX 


This  subroutine  permits  the  calculation  of  the  refractive  index,  its 
spatial  derivatives,  and  the  absorption  coefficient,  all  as  functions 
of  the  local  atmosphere  and  its  state  of  ionization.  This  subroutine 
also  permits  the  output  of  "debugging  data"  which  is  called  the  R 
vector  whose  components  are  defined  in  Table  3.  See  subroutine 
OUTPUT  for  additional  information. 
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R(D  =  X 

R<2)  =  PXR  =  |5. 

R(3)  =  PXTHET  = 

R(4)  =  PXPHI  =  -g- 

R(5)  =  Y 

R(6)  =  YR  =  Yf 

R(7)  =  YTHETA  =  YQ 

D 

R(8)  =  YPHI  =  Y? 

R(9)  =  PYR  =  8Y/8r 

R(  1 0)  =  PYTHET  =  8Y/80 
R(1 1)  =  PYPHI  =  8Y/8«R 
R(12)  =  Z 

R(13)  =  PZR  =  8  Z/8r 
R(  14)  =  PZTHET  =  8Z/86 
R(1 5)  =  PZPHI  =  QZ/8q> 

R(16)  =  COSPSI  =  cos  if 
R(17)  =  SINPSI  =  sin  t 
R(18)  =  YSI  =  Y  sin  + 

R(19)  =  YCI  =  Y  cos  I 
R(20)  =  TE  1  =  (1  -  X)2  -  Z2 

Table  3.  Nomenclature  Describing  the  R  Vector. 
(Page  1  of  4) 
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R(21 )  =  TE2  =  (YSI)4  +  4TE1  (YCI)2 
R(22)  =  TE3  =  8(YCI)2  Z(1  -  X) 

R(23)  =  R2S  =  Rg2 

R(24)  =  R  IS  =  Rg 

R(25)  =  THET2S  =  20g 

R(26)  =  THETls  =  0g 

R(27}  s  SI  s  Sj 

R(28)  =  S2  =  S2 

R(29)  =  D1  =  dl 

R(30)  =  D2  =  d2 

R(31)  =  TE4  =  d*  +  d22 

R(32)  =  TE5  =  2x[zd1  +  (1  -  X)d2]/TE4 

R(33)  =  TE6  =  1  -  [ 2X(1  -  X)d{  -  Zd£  ]/TE4 

R{34)  =  R2M  =  Rm2 

R(35)  =  RIM  =  Rm 

R(36)  =  THET2M  =  20M 

R(37 )  =  THET1M  =  ©M 

R(38)  =  AMI  =  Mj 

R(39)  =  AM  2  =  M2 

R(40)  =  TE7  =  MJd1  -  M^ 

Table  3.  Nomenclature  Describing  the  R  Vector. 

(Rage  2  of  4) 
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R(4 1 )  =  TE8  =  M1d2  +  M2d1 
R(42)  =  AO  =  aQ 

R(43)  =  BO  =  b 

o 

R(44)  =  TE9  =  Sl2  +  S22 
R(45)  =  A4  =  a4 

R(46)  =  B4  =  b4 

R(47 )  =  A5  a5 
R(48)  =  B5  =  b& 

R(49)  =  PNPX  =  8n/8x 

R(50)  =  TE 10  =  (sint)2(YSI)2  +  2TE1  cos2i|r 

R(51)  =  TE11  =  4Z(1  -  X)cos2i 

R(52)  =  A6  =  a6 

R(53)  =  B6  =  b6 

R(54)  =  PNPY  =  8»i/  8Y 

R(55)  =  A7  =  a? 

R(56)  =  B7  =  b? 

R(57)  =  PNPZ  =  8n/8Z 
R(58)  =  A8  =  ag 
R(59)  =  B8  =  bg 

2  1  2 

R(60)  =  TE  1 3  =  1/(2  co»  0  ♦  J  «in  6) 

Table  3.  Nomenclature  Describing  the  R  Vector. 

(Page  3  of  4) 
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R(61)  =  A1  =»  a 

R(62)  =  B1  »  b 

R(63)  =  A2  =  a2 

R(64)  =  B2  =  b2 

R(65)  =  TEH  =  (W^2  YSI 

R(66)  =  PPSIPT  =  at/89 

R(67)  =  PPSIPR  =  3t/8r 

R(68)  =  PPSIPP  =  at/aq> 

R(69)  =  TE 14  =  J~arZ  +  cQ2  +  a* 


Table  3.  Nomenclature  Describing  the  R  Vector. 
(Page  4  of  4) 
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I.  SUBROUTINE  ELECTX 


The  whole  process  of  determining  the  free  normalised  electron 
density  and  its  spatial  derivatives  at  any  point  in  space  is  explored 
in  this  subroutine.  For  the  case  under  consideration  this  is  a 
simple  sphere  where  the  values  in  Equation  74  are  represented  by 
the  following  components  of  the  W  vector 

W(28)  =  &  W (34)  »  A;  W(35)  »  n 
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J.  SUBROUTINE  BIGR 


This  subroutine  is  used  to  calculate  the  distance  from  the  center  of 
the  ionising  source  to  the  spatial  point  r,  6,  V  at  which  the  electron 
density  and  spatial  gradients  are  desired.  For  the  simple  sphere 
this  distance  is  W(28)  *  In  addition  to  this  it  calculates  the 
angle  that  l€  makes  with  respect  to  the  vertical  distance  passing 
through  the  center  of  the  ionising  source. 
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K.  SUBROUTINE  MAGY 


The  normalised  earth's  magnetic  field  and  its  spatial  derivatives 
are  calculated  in  this  subroutine.  In  this  program  it  is  represented 
by  an  assumed  magnetic  dipole  field  with  the  dipole  located  at  the 
center  of  the  earth. 
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L.  SUMOUTINC  C0LFR2 


The  normalised  collision  frequency  and  its  derivatives  at  a  spatial 
point  are  determined  in  this  subroutine.  In  this  presentation  it  is 
assumed  that  collision  frequency  varies  exponentially  with  height. 
One  should  refer  to  the  definition  of  the  W  vector  where  the  required 
coefficients  for  the  height  stratifications  are  defined. 
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M.  SUBROUTINE  RCOORD 


This  subroutine  permits  the  transformation  from  the  earth  centered 
geomagnetic  spherical  coordinate  system,  to  the  radar  coordinate 
system. 
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N.  SUBROUTINE  OUTONE 


In  this  presentation  this  subroutine  is  used  to  list  the  input  data  and 
the  initial  conditions  that  define  the  ray  tracing  problem  on  Tape 
Unit  6.  Table  5-1  illustrates  this  output  by  the  statement  of  the  in¬ 
put  RECORD  and  the  next  ninety-one  words  of  data  which  define  in 
order  the  first  70  components  of  the  W  vector,  followed  by  the  first 
21  components  of  the  V  vector.  This  subroutine  can  be  greatly  im¬ 
proved. 
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O.  SUBROUTINE  OUTPUT 


This  subroutine  will  presently  yield  the  output  data  on  Tape  Unit  6 
that  is  summarized  in  Tables  5  and  6.  Following  the  information 
that  is  the  result  of  subroutine  OUTONE,  Tables  5-1  and  5-2 
illustrate  the  information  obtainable  after  each  integration.  Begin¬ 
ning  with  the  first  word  this  information  has  the  following  meaning: 

Integration  number  W(70) 

Vector  components  V(2),  V(3) 

Height  above  earth  surface  (km) 

Angle  8  in  degrees 
Angle  9  in  degrees 

Vector  components  V(7)  through  V(21) 

Vector  components  W(l)  through  W(ll) 

Vector  components  W(28),  W(29).  W(38).  W(39) 

Vector  components  W(55)  through  W(70) 

Vector  components  W(80)  through  W(90) 

Under  certain  conditions  of  vector  components  W(68)  and  W(69)  the 
R  vector  described  in  subroutine  RINDEX  will  be  listed  between 
each  of  these  sets  of  data  for  each  time  that  the  RINDEX  subroutine 
is  entered. 

In  addition  to  this  data,  the  data  summarized  in  Table  6  is  listed  on 
Tape  Unit  10. 
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P.  INPUT-OUTPUT 


Table  4  illustrates  the  format  that  is  necessary  for  the  input  data. 
Tables  5  and  6  illustrate  the  format  of  the  output  data  that  is  cur¬ 
rently  obtained  by  use  of  Subroutines  OUTONE  and  OUTPUT. 
These  two  subroutines  can  be  easily  modified.  As  they  are  pre¬ 
sented  in  this  report  their  primary  purpose  was  to  collect  "debug¬ 
ging"  data. 


199 


i 

£ 


i 


D 

D 

D 

D 

0 

D 

D 

n 

n 

r 


i 

i 

i 

[ 


i  ia  Table  4.  Inout  Data  for  a  Sohorical  tononharo. 


RM  61TMP-32 


Tabl*  5-1 .  Ouhaif  Data  far  "Dabuoaina"  Rtnow  from  Tom  6. 
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